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ABSTRACT 

A  radar  target,  acting  as  a  scatterer  of  an  incident 
electromagnetic  wave,  can  be  considered  as  a  linear  time- 
invariant  system.  Previous  work  has  shown  that  the  target's 
pole  locations  are  independent  of  the  incident  electromagnetic 
excitation,  including  incident  wave  shape,  aspect  and 
polarization.  This  thesis  develops  the  Kumaresan-Tuf ts  and 
Cadzow-Solomon  signal  processing  algorithms  into  computer 
routines  and  evaluates  their  pole  extraction  performance. 
Data  used  to  evaluate  the  extraction  algorithms  includes 
synthetic  and  integral  equation  generated  signals  with 
additive  noise,  in  addition  to  measurements  of  scattering  by 
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I.  INTRODUCTION 

A  radar  target,  acting  as  a  scatterer  of  a  specified 
incident  electromagnetic  wave,  can  be  considered  as  a  single 
input,  single  output,  linear  time-invariant  (LTD  system  for 
a  fixed  field  observation  point.  The  target  can  thus  be 
considered  as  a  transfer  function  with  poles  and  zeros.  Baum 
demonstrated  at  the  Air  Force  Weapons  Laboratory  that  a 
target's  induced  current  response  to  an  incident  electro- 
magnetic wave  has  identifiable  poles  determined  by  the 
composition  and  structural  geometry  of  the  target  [1]  .  In 
1974,  Moffatt  and  Mains  proposed  that  the  target's  scattered 
field  pole  locations  are  independent  of  the  incident 
electromagnetic  excitation,  including  aspect  and  polarization 
[2].  Morgan  has  proven  theoretically  that,  for  the  case  of 
a  conducting  target,  the  scattering  response  contains  complex 
natural  resonances  which  are  independent  of  the  incident 
electromagnetic  excitation  [3] .  By  determining  the  poles  of 
a  target's  response,  aspect  independent  target  identification 
can  be  accomplished  through  the  use  of  electromagnetic  natural 
resonances . 

Although  the  concept  of  radar  target  identification 
through  the  use  of  natural  resonances  was  first  proposed  in 
1974  by  Mains  and  Moffatt   [2]  ,   only  recently  have  signal 


processing  techniques  been  applied  to  locate  the  poles  in  a 
radar  target's  response  in  the  presence  of  noise.  Kumaresan- 
Tufts  [4]  and  Cadzow-Solomon  [5]  have  each  developed 
algorithms  which  have  proven  successful  in  the  presence  of 
noise.  This  thesis  develops  computer  routines  based  upon 
these  two  algorithms  and  examines  their  respective  performance 
and  appropriateness  using  a  variety  of  scattering  data. 

A.   THE  PROBLEM 

Since  the  performance  of  signal  processing  methods  varies 
under  different  conditions,  a  system  employed  to  identify 
targets  would  possibly  reach  a  decision  based  on  the  combined 
output  of  several  signal  processing  methods.  For  example,  the 
Kumaresan-Tuf ts  and  Cadzow-Solomon  methods  could  be  used  to 
extract  poles  from  the  response  of  scale  model  targets.  The 
information  so  gathered  could  be  used  to  build  a  data  base  for 
comparison  with  data  similarly  obtained  in  actual  field  use. 
The  results  of  this  system  would  serve  as  one  input  to  a 
larger  system.  Other  methods  would  provide  input  to  the 
system,  such  as  the  K-pulse  method  of  Kennaugh  [6]  and  the 
annihilation  filter  used  by  Dunavin  [7] ,  Morgan  and  Dunavin 
[8]  and  Chen  [9],  As  the  name  suggests,  an  annihilation 
filter  annihilates  the  target's  poles.  A  system  using  the 
annihilation  filter  concept  would  contain  many  such  filters, 
each  previously  designed  to  cancel  the  poles  of  a  specific 


known  target.  In  actual  field  use,  a  radar  target's  response 
would  be  input  into  each  of  the  filters,  and  the  target 
selected  would  be  that  matching  the  filter  whose  output 
exhibits  the  lowest  signal  energy. 

A  system  used  to  identify  radar  targets  would  require  the 
following  concept  of  employment.  First,  information  required 
by  each  of  the  sub-systems  would  be  obtained  for  every  target 
class  of  concern.  In  actual  field  use,  this  information  would 
be  compared  against  actual  radar  target  responses.  The  system 
would  then  determine  the  identity  of  the  target  based  on  the 
input  from  each  of  its  sub-systems. 

B.   BACKGROUND 

Consider  a  perfectly  conducting  target  illuminated  by  an 
electromagnetic  field.  The  current  induced  on  the  surface  of 
this  target  at  a  given  point  must  satisfy  the  magnetic  field 
integral  equation  (MFIE) , [10] 

J(r,0=2nxH'(r,0+JJ  K(r,r  ,t) J(r/~1  ?c"r  l)dS         (1) 

where  n  is  an  outward  unit  vector  normal  to  the  surface  of 
the  object,  J  is  the  surface  current  density,  .H  .  is  the 
incident  magnetic  field,  and  K  is  a  Green's  function  dyadic. 
The  entire  equation  is  most  easily  understood  as  the  sum  of 
driven  currents  and  "feedback"  currents  corresponding  to  the 


cross-product  term  and  surface  integral  term  respectively. 
The  term  driven  by  the  magnetic  f  ield,  2n*H  i  forms  the  physical 
optics  portion  of  the  total  current.  Physical  optics 
describes  the  cross-product  term  as  the  induced  current 
without  interaction  with  the  rest  of  the  body.  The  Green's 
function  kernel  describes  the  current  at  a  point  on  the  object 
due  to  the  feedback  of  currents  from  every  other  point  on  the 
object,  as  previously  illuminated  by  the  incident  field.  The 
current  at  each  point  is  then  summed  over  the  surface  of  the 
object.  Note  that  the  surface  integral  term  is  of  principal- 
value  type;  the  integral  excludes  the  point  r=r  . 

Once  the  incident  magnetic  field  is  no  longer  present, 
the  solutions  of  (1)  are  considered  the  natural  modes  of  the 
object.  These  natural  modes  are  of  the  form,  Jnexp(sn)  .  The 
natural  resonance  frequencies   sn   are  of  the  form, 

sn=(V:K  (2) 

where  on  is  the  damping  rate  in  Nepers/sec  and  wn  is  the 
frequency  in  radians/sec.  The  natural  resonances  of  (2)  are 
functions  of  the  structural  geometry  of  the  object  and  are 
independent  of  the  incident  magnetic  field.  To  understand 
how  these  natural  resonances  are  unique  to  the  geometry  and 
composition  of  the  object,  consider  a  set  of  points  on  the 
object  previously  illuminated  by  the  incident  field,  so  that 
H  =0 .   The  current  at  a  given  point  in  the  set  is  due  to  the 


infinite  number  of  feedback  currents  from  every  other  point 
in  the  set.  Recall  that  these  feedbacks  are  described  by  the 
Green's  function  kernel  in  the  integral  term  of  (1).  Since 
the  set  of  points  previously  illuminated  is  physically  located 
on  the  same  object,  the  infinite  number  of  paths  that  connect 
a  point  with  all  other  points  in  the  set  is  the  same  for  all 
points  in  the  set.  The  infinite  number  of  paths  are  unique 
to  the  structural  geometry  of  the  object  and  correspond 
exactly  to  the  infinite  number  of  paths  taken  by  currents 
which  feedback  to  a  given  point  via  the  Green's  function 
kernel.  Finally,  the  composition  of  the  target  determines  the 
surface  current  density  on  the  object.  Although  an  infinite 
number  of  resonances  exists  in  any  object,  only  a  limited 
number  of  these  will  be  measurably  excited  by  an  incident 
field  of  finite  bandwidth.  These  resonances  described  in  (2) 
appear  as  complex  conjugate  pairs  in  the  left-half  portion  of 
the  s-plane. 

In  the  far-field,   the  back-scattered  response  of  a  target 
to  an  incident  plane  wave  is  of  the  form 


S^-rP't)=4^F  a|jJjP*J(r,t-f-r7c)dS' 


(3) 


where  c  is  the  speed  of  light  and  p  is  the  unit  vector  whose 
direction  matches  that  of  the  plane  wave's  propagation. 


Equation  (3)  is  the  result  of  integrating  the  current  at 
each  point  on  the  target  surface  for  a  fixed  point  in  the  far- 
field.  Recall  that  the  current  at  each  point  on  the  target 
is  defined  by  (1).  Thus,  the  back-scattered  far-field  can 
be  obtained  by  substituting  (1)  into  (3) : 


H(-rp,t)=u(t-r/C)fHpJ-rp,t)+J_Jln(-rp,t)exp(snt)  \        (4) 


The  currents  in  (1)  produce  the  field  in  (4).  In  fact,  each 
term  in  (4)  corresponds  to  the  term  in  (1)  which  produced  it. 
Specifically,  the  first  term  in  (4)  describes  the  physical 
optics  scattered  field  generated  by  the  2nxH  current  which, 
of  course,  is  the  first  term  in  (1).  Similarly,  the  second 
term  in  (4)  is  produced  by  the  source-free  currents  defined 
by  the  second  term  in  (1) .  Like  the  current  described  in  (1) , 
the  field  in  (4)  is  the  sum  of  two  terms,  a  driven  term  and 
a  term  containing  feedbacks. 

The  results  of  (4)  can  also  be  seen  as  two  forms  of  the 
Singularity  Expansion  Method  (SEM)  developed  by  Baum  [1] .  As 
shown  by  Morgan  [10] ,  during  the  early-time  portion  of  the 
target's  response,  the  scattered  field  is  composed  of  the 
physical  optics  scattered  field  and  a  "Class  2"  form  of  the 
SEM  expansion.  The  class  2  SEM  expansion  corresponds  to  the 
second  term  of  (4),  wherein  the  coefficients  Hn  are  time- 
varying  as  the  wave  passes  over  the  target,  since  the  currents 


producing  this  portion  of  the  field  are  integrated  over  a 
time-varying  surface  area.  At  the  instant  the  wave  passes  the 
last  point  of  the  target,  the  physical  optics  field  vanishes 
and   the   remaining  term  in   (4)   is  produced  by  constant 

coefficients  H  .   The  coefficients  H   are  constant  at  this 

11 .  n 

instant  since  the  surface  area  in  the  integral  in  (3)  is  now 
constant.  This  instant  also  marks  the  transition  of  (4)  from 
a  "class  2"  SEM  expansion  of  time-varying  coefficients  to  a 
"class  1"  SEM  expansion  of  constant  coefficients.  The 
scattered  field  due  to  a  plane  wave  is  therefore  composed  of 
a  physical  optics  term  and  a  class  2  SEM  expansion  in  the 
early-time,  and  a  simple  class  1  expansion  in  the  late-time. 

Actual  measurement  of  the  scattered  far-zone  field  would 
be  greatly  aided  by  knowledge  of  the  transition  time  of  the 
field  from  early  time  to  late  time.  From  [10] ,  this 
transition  for  a  monostatic  radar  would  occur  at  At=T+2 (D+d)/c 
seconds  after  radar  turn-on.  Here,  T  is  the  pulse  duration, 
D  is  the  target's  dimension  along  the  direction  of  wave 
propagation,  d  is  the  distance  between  the  target  and  the 
measurement  point  and  c  is  the  speed  of  light. 

The  discussion  presented  in  this  section  was  extracted 
from  work  done  by  Morgan  in  [10] .  The  reader  is  referred  to 
this  work  for  a  more  detailed  treatment  of  the  material  in 
this  section. 


C.   HISTORY 

The  results  of  the  previous  section  form  the  basis  for  the 
hypothesis  that  the  natural  resonances  found  in  the  scattering 
response  of  a  target  to  an  incident  electromagnetic  wave  are 
unique  to  that  target.  Additionally,  only  a  finite  set  of 
these  natural  resonances  are  measurably  excited  by  a  wave  of 
finite  bandwidth.  In  1974,  Moffatt  and  Mains  proposed  that 
the  extraction  of  resonances  from  a  target's  response  to 
electromagnetic  excitation  could  be  used  for  target 
identification.  This  work  related  to  earlier  work  in  1965, 
when  Kennaugh  and  Moffatt  first  developed  the  concept  of  a 
radar  target  as  a  linear  time  invariant  system.  Poles  in  the 
z-plane  are  directly  related  to  the  natural  resonances  of  a 
target 

z  =e 

(5) 

where  sn  is  given  by  (2)  and  At  is  the  sampling  interval  in 
seconds.  Hence,  pole  extraction  involves  resonance 
identification.  The  use  of  pole  extraction  algorithms  is 
discussed  in  the  next  chapter. 


II.   POLE  EXTRACTION  ALGORITHMS 

The  use  of  pole  extraction  algorithms  to  identify  radar 
targets  is  discussed  in  this  chapter.  A  brief  discussion  of 
two  methods  precedes  the  in-depth  evaluation  of  the  Kumaresan- 
Tufts  and  Cadzow-Solomon  algorithms.  The  evaluation  of  the 
latter  two  algorithms  occurs  in  two  stages.  First,  each 
algorithm  will  be  evaluated  in  its  ability  to  extract  poles 
from  data  with  known  poles.  Some  of  the  data  processed  was 
generated  at  various  signal  to  noise  ratios  by  a  computer 
program  written  by  Morgan  [11] .  Additional  data  was  produced 
by  Morgan's  time-domain  thin  wire  integral  equation  computer 
program  [12].  In  the  second  stage,  a  side  by  side  comparison 
is  made  of  poles  extracted  by  each  method  using  transient 
scattering  measurements  for  a  thin  wire  and  for  various  model 
aircraft.  Comparisons  between  the  two  methods  are  made  as  the 
aspect  of  the  aircraft  is  varied. 

A.   PREVIOUS  WORK 

1.   Direct  Minimization 

The  most  direct  way  to  determine  the  natural 
resonances  in  a  target's  response  is  to  minimize  the  mean- 
square  error  between  the  modeled  signal  and  the  received 
signal.   In  [10],  Morgan  determined  that  the  late-time  target 


response  to  a  radar  could  be  represented  as  a  sum  of  damped 
sinusoids  given  by 


00    o   t 
y(t)  =  lA1e  '  cos(oj,t+e,)  (6) 

i  =  i 


The  frequency,  o^  ,  and  damping  rate,  0  ,  are  the  same 
parameters  found  in  the  natural  resonance  defined  in  (2) . 
Phase,  0  ,  and  amplitude,  A,,  are  the  remaining  parameters. 
The  representation  in  (6)  is  the  sum  of  an  infinite  number  of 
resonances.  The  sampled  response  to  an  incident  wave  of 
finite  bandwidth  can  be  modeled  as 


y(nat)  =  y  =  £  A,e  l"  cos  (o^nat+G,)  (7) 

n  1  =  1 


where  At  is  the  sampling  interval  in  seconds.  The  four 
parameters  of  (7)  must  be  adjusted  to  minimize  the  sampled 
mean-square  error  signal 

eHyn-yn>2  (8) 

between  the  actual  discrete  sampled  received  signal  y  and  the 

n 

modeled  signal  y  .  The  processing  required  in  this 
minimization  problem  is  both  inefficient  and  highly  non- 
linear. Nevertheless,  Chong  used  this  method  to  process 
mathematically-generated  data  down  to  15.0  dB  signal-to-noise 
(SNR)  ratio  [13]  . 
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2.   Prony's  Method 

As  in  direct  minimization,  Prony's  approach  to 
resonance  classification  focuses  on  the  late-time  portion  of 
a  radar  target's  response.  However,  linear  processing  and 
root  solving  are  used.  The  late-time  response  is  modeled  as 
the  output  of  an  LTI  system  of  order  KD.«  Each  signal  received 
at  some  discrete  sample,  n,  is  considered  to  be  the  weighted 
sum  of  KD  previous  signals.  Thus,  the  finite  term 
approximation  of  the  received  late-time  signal,  y  ,  is  defined 

n 

by 


y  =lhiY 


i=i 


(9) 


The  z-transform  of  (9)  is 

K„-i  . 


2^-b1ZK^1-b2ZK^1...-bKD=0 


(10) 


The  roots  of  this  polynomial  in  z  are  the  poles  of  the  system 
model.  Therefore,  the  key  to  extracting  the  poles  in  the 
system's  response  lies  in  solving  for  the  coefficients  b,  of 
(9)  . 

A  set  of  KD+M  received  signals  in  M  equations  (9)  can 
be  arranged  in  matrix  form  as 


KD-1 


v.---y 


KD+M-2 


rbv 


y„     n 


_  ■*  KD+M-1 


(11) 
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In  Prony's  original  method,  the  data  matrix  d  is 
exactly  determined,  and  the  coefficient  vector,  b,  is  solved 
using  linear  computations.  In  the  presence  of  noise,  Prony 
overdetermines  the  data  matrix  by  setting  M>KD  and  solves  for 
the  coefficient  vector  by  obtaining  the  least-squares  solution 
to  the  system  of  equations. 

The  Prony  method  has  two  major  problems.  First,  the 
poles  obtained  by  the  least  squares  solution  to  the 
overdetermined  matrix  may  be  strongly  perturbed  by  noise  [14]  , 
since  noise  does  not  satisfy  the  causal  model  of  the  system. 
Second,  the  order  of  the  system  is  generally  not  known  a 
priori.  When  the  estimated  order  is  greater  than  the  actual 
order,  poles  due  to  noise  are  generated.  Prony's  method 
offers  no  technique  for  distinguishing  between  the  signal 
poles  and  the  extra  poles  caused  by  overestimation  of  the 
system's  order.  If  the  estimated  system  order  is  less  than 
the  actual  order,  actual  poles  are  lost  and  the  remaining 
poles  are  perturbed  from  their  true  positions. 

B.   KUMARESAN-TUFTS  ALGORITHM 

The  Kumaresan  and  Tufts  pole  extraction  algorithm  was 
developed  by  adapting  Prony's  method  to  reduce  the  problems 
addressed  in  the  preceding  section.  The  Kumaresan-Tuf ts 
algorithm  modifies  the  least-squares  Prony  method  in  three 
ways : 
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1.  Processed  signals  are  arranged  in  a  data  matrix 
based  on  a  non-casual  model  of  the  system. 

2.  The  model  of  the  system  is  deliberately 
overestimated. 

3.  The  system  of  equations  determined  by  the  above 
two  criteria  is  solved  by  using  singular  value 
decomposition  (SVD) . 


Kumaresan  demonstrates  in  [15]  that  the  use  of  singular 
value  decomposition  tends  to  force  the  extra  poles  of  the 
excess-order  system  inside  the  unit  circle,  while  the  non- 
causal  arrangement  of  the  signals  tends  to  force  the  signal 
poles  outside  the  unit  circle.  The  excess  order  of  the  system 
model  reduces  the  effects  of  noise  on  the  actual  poles.  Since 
the  noise  is  stationary  and  stable,  it  looks  the  same  in 
forward  and  backward  time. 

1 .   Equations 

Recall  that  in  (9),  Prony ' s  technique  defines  the 
received  late-time  signal  as  the  weighted  sum  of  kd  previous 
signals,  where  KD  is  presumed  to  be  the  order  of  the  system. 
Kumaresan  models  the  same  late-time  signal  as  the  weighted  sum 
of  KD  future  signals,  where  KD  is  greater  than  the  estimated 
order  of  the  system.   This  non-casual  model  is  given  by 


KD 


y=Ib',y  (12) 
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A  system  of  M  such  prediction  equations  can  be  written  in 
matrix  form  as 


r  y 


KD+M-1 


Or,  in  matrix  notation, 


\a\  -[l] 


(13) 


Dy-b=y 


(14) 


As  in  Prony '  s  method,  the  coefficients  b',  are  coefficients  of 
a  polynomial  in  z  that  models  the  system's  late-time  response. 
Two  simple  manipulations  of  either  data  matrix  leads  to  the 
relationship  between  the  coefficients  of  the  Prony  model  and 
the  prediction  coefficients  of  the  Kumaresan-Tuf ts  model. 
With  b0=-l  ,  a  prediction  coefficient  is  related  to  an 
autoregressive  coefficient  by 


w  -  i-i 


(15. 


From  the  above  relationship,  it  can  be  shown  that  the  complex 
pole  pairs  of  the  causal  model  are  merely  conjugate 
reflections  across  the  unit  circle  of  the  pole  pairs  in  the 
non-causal  model. 

2.   Singular  Value  Decomposition 

The  non-causal  arrangement  of  late-time  signals  in  a 
set  of  system  equations,  and  subsequent  processing  through 
singular  value  decomposition,  combine  to  separate  the  signal 
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and  noise  into  orthogonal  spaces.  As  discussed  in  the 
preceding  paragraph,  poles  of  the  non-causal  model  are 
reflected  outside  the  unit  circle.  Kumaresan  demonstrates  in 
[15]  that  the  extra  poles  of  the  excess-order  system  can  be 
forced  inside  the  unit  circle  through  the  use  of  SVD . 

Singular  value  decomposition  factors  the  i  MXKD  data 
matrix  D  into  the  product  of  the  matrices: 


D  =USVT 


(16) 


The  columns  of  U  (MXM)  are  eigenvectors  of  DyDy  an<^  tne 
columns  of  V  (KDXKD)  are  eigenvectors  of  yDy  •  If  r  is  the 
rank  of  the  data  matrix,  D  ,  the  diagonal  matrix  I  (MXKD) 
contains  r  singular  values  which  are  the  square  roots  of  the 
nonzero  eigenvalues  of  both  D^D  and  DyDy.-  Bv  rearranging 
the  three  matrices  in  the  product,  the  pseudoinverse  of  Dy  can 
be  obtained  as 

D!=VI+UT 
y  (17) 

where  i+  is  a  (KDXM)  matrix  whose  singular  values  on  the 
diagonal  are  the  reciprocals  of  those  in  the  I  matrix. 
Finally,  the  coefficient  vector  b+,  of  minimum  Euclidian  norm, 
is  given  by 

b-=D;y  <18> 
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The  coefficient  vector  fc>+  so  obtained  is  the  minimum  length 
least-squares  solution  to  (14).  In  other  words,  b+  is  the 
best  possible  solution  to  (14)  .  In  the  case  of  noiseless 
data,  the  extraneous  poles  generated  by  the  excess-order  model 
will  always  be  inside  the  unit  circle  when  b+  is  used.  This 
result  is  generally  true  for  noisy  data. 

3.  Bias  Compensation 

Kumaresan  and  Tufts  [4]  observed  that  the  addition  of 
noise  perturbed  the  singular  values  of  the  z  matrix  of  (16). 
If  the  perturbation  of  these  singular  values  is  not 
compensated,  both  the  signal  poles  and  extraneous  poles  are 
biased  towards  the  unit  circle.  Kumaresan  and  Tufts  used  a 
compensation  method  which  reduced  the  bias  in  their  work,  but 
did  not  derive  an  analytical  justification.  In  [16],  Norton 
derived  a  more  valid  bias  compensation  method  based  on  the 
eigenvalue  shifting  theorem. 

4.  Kumaresan  and  Tufts  Compensation 

If  the  actual  order  of  the  system  is  K^  _,  then  the 
first  Kp  singular  values  of  the  I  matrix  in  (16)  are  non- 
zero. The  remaining  KD-Kp  singular  values  are  considered 
noise  singular  values  and  are  zero  in  the  case  of  noiseless 
data.  The  addition  of  noise  perturbs  the  first  KD  signal 
singular  values  and  increases  the  noise  to  some  non-zero 
value.  Kumaresan  and  Tufts  compensated  for  this  increase  in 
the  singular  values  due  to  the  noise  by  subtracting  the 
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average  of  the  noise  singular  values  from  the  signal  singular 

values.   The  noise  singular  values  were  then  set  to  zero. 

5.   Compensation  Based  on  Eigenvalue  Shifting  Theorem 

As  described  in  the  previous  section,  the  singular 

values   of   the  matrix   Dy   are   the   square  roots   of   the 

eigenvalues  of   D„Dy,  and  dtD  •   Assume  the  noisy  data  matrix 

y  y-       y  y  ' 

can  be  represented  by  d  =s+N'»  where  N  is  composed  of  the  wide- 
sense  stationary  white  noise  process  v.,  given  by 


Ny= 


Vv  .  .  .  v  ' 


M+K 


D  _l 


(19) 


The  expected  value  of  D  D^  can  be  obtained  by 


y  y 


DyDy=E[  (S+N)  (S+N)T]=E[SST]+E[SNT]+E[NST]+E[NNT]      (20) 

Since  S  is  deterministic,  E[SST]=SST.  Assuming  the  noise  is 
zero  mean,  the  two  cross  products  are  zero.  Because  we  assume 
the  noise  is  wide-sense  stationary  and  white,  E[NNT]=a^I, 
where  a2  is  the  noise  variance  and  I  is  the  identity  matrix. 
The  expected  value  of  ]  D  DT   thus  becomes 

y   y- 


E[DyDTy]=SST+o2vI 


(21) 
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Similarly,  the  expected  value  of  DTD  ,     the  other  source  of 

y  y 

singular  values,  is 

E[DTyDy]=STS+o2vI  (22) 

The  assumption  in  the  results  of  (21)  and  (22)  is  that  the 
diagonals  of  E[N  N]=E[NN  ]  equals  the  noise  variance  <j2  .• 
Equations  (21)  and  (22)  show  that  in  the  mean,  the  squares  of 
the  singular  values  of  D  are  increased  by  the  noise  variance. 
The  results  lead  to  the  method  of  eigenvalue 
compensation  recommended  by  Norton  in  [16].  Recall  from  (16) 
that  the  eigenvalues  of  D  are  on  the  diagonal  of  the  I 
matrix  returned  by  the  singular  value  decomposition  of  D  • 
If  KD'  is  the  actual  order  of  the  system,  and  kd  is  the 
estimated  order  of  the  system  then  the  remaining  KD-KD' 
singular  values  of  the  I   matrix  can  be  squared  and  averaged 

2 

to  obtain  an  estimate  of  the  noise  variance,  °v  .  These  noise 
singular  values  can  then  be  set  to  zero.  The  first  KD 
singular  values  of  the  I  matrix  are  then  squared  and  reduced 
by  subtracting  the  estimate  of  the  noise  variance.  The  square 
root  of  the  difference  becomes  the  new  first  KD  singular 
values  of  the  compensated  Z  matrix.  Calculations  according 
to  (17)  and  (18)  can  then  be  carried  out  in  a  normal  manner 
to  obtain  poles  in  the  presence  of  the  noise.  Eigenvalue 
compensation  requires  an  estimate  of  the  actual  order  of  the 
system.  Methods  to  obtain  this  estimate  are  discussed  in 
Chapter  III. 


6.   Performance 

The  Kumaresan-Tufts  algorithm  was  programmed  in 
Fortran  and  tested  on  various  types  of  data.  The  program 
appears  in  Appendix  A. 

a.   Synthetically  Generated  Data 

The  starting  point  for  evaluating  the  performance 
of  the  Kumaresan-Tufts  algorithm  was  with  synthetically 
generated  data  of  the  form  given  by  (8)  and  shown  here  again 
for  convenience 


y=XA,e  '   cosfw.nAt+G,)  (8> 

n  1=1 


Again,   A^o^w,^  ,  are  the  amplitude,  damping  rate,  frequency 
and  phase  of  a  set  of  N  damped  sinusoids.    Noisy  data  was 
created  by  adding  stationary  white  noise. 
1.   Noise  Performance 

The  algorithm  was  evaluated  at  various  SNR's, 
ranging  from  90.0  dB  to  7.0  dB .  These  SNR's  are  ratios  of 
signal  energy  to  noise  energy  rather  than  the  ratio  of  signal- 
to-noise  power.  Synthetic  data  so  generated  more  closely 
resembles  the  exponential  decay  of  signal  power  typical  in 
actual  radar  measurements. 
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Figure  1  shows  the  signal  produced  by  two  s- 
plane  poles  at  90.0  dB.  Figures  2  through  6  depict  the  poles 
extracted  from  this  signal  at  SNR ' s  ranging  from  90.0  dB  to 
7.0  dB .  Obtained  poles  are  shown  at  their  positions  within 
the  upper  right  hand  quadrant  of  the  unit  circle  in  the  z- 
plane.  Not  shown  are  conjugates  of  each  pole  which  are 
located  below  the  real  axis  outside  the  figure  boundaries. 

Figures  2  through  6  demonstrate  outstanding 
performance  on  noisy  data,  even  at  SNR '  s  of  7.0  dB .  The 
scaling  needed  to  show  a  discernible  difference  between 
results  obtained  at  30.0  dB  and  7.0  dB  would  necessarily 
exclude  one  of  the  poles  from  the  enlarged  figure.  The 
average  distance  of  the  trial  poles  obtained  in  the  7.0  dB 
SNR  signal  from  the  true  poles  is  on  the  order  of  10"3-  This 
magnitude  corresponds  to  that  of  the  average  estimate  of  the 
noise  variance  obtained  in  successive  trials  with  this  signal. 
The  correlation  between  the  distance  of  trial  poles  from  true 
poles  and  the  noise  variance  estimate  was  consistently 
observed  with  each  of  the  different  signal-to-noise  ratios 
used.  Figure  7  depicts  the  signal  of  Figure  1  severely 
corrupted  by  noise  having  7.0  dB  SNR. 

As  discussed  previously,  the  signal-to-noise 
ratio  used  in  the  synthetically  generated  data  is  the  ratio 
of  energy.  Figure  8  depicts  the  results  of  pole  extraction 
from  the  signal  shown  in  Figure  7,  but  with  a  late-time 
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Figure   1.   Signal  Containing  two  S-Plane  Poles, 


9  0.0  dB  SNR 
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Figure  3.    Kumaresan-Tufts  Poles,  Synthetic  Data,  30.0dBSNR 
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Figure  6.   Kumaresan-Tufts  Poles,  Synthetic  Data,  7 .  0  dB  SNR 
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Figure   8.   Kumaresan-Tuf ts  Pole  Extraction,  7.0  dB  SNR 
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beginning  ten  nanoseconds  later.  Since  the  SNR  is  calculated 
over  twenty  nanoseconds  for  both  signals,  the  signal  power  at 
some  later  time  will  clearly  be  less  than  the  power  ten 
nanoseconds  earlier.  The  results  in  Figure  8  show  complete 
breakdown  of  the  algorithm's  ability  to  extract  poles.  The 
trial  poles  shown  are  the  poles  closest  to  the  true  poles,  and 
yet  they  are  located  at  positions  whose  reflections  are  inside 
the  unit  circle  where  noise  poles  are  typically  located. 

The  preceding  results  show  outstanding 
accuracy  for  full-length  noisy  data  but  a  complete  breakdown 
of  the  algorithm  for  the  same  signal  with  a  later  transition 
to  late-time.  These  initial  observations  are  supported  by 
similar  findings  presented  in  this  thesis. 

b.   Thin  Wire  Integral  Equation  Generated  Data 

For  simple  objects  such  as  a  thin  wire,  the  radar 
response  of  that  object  can  be  computed  by  establishing 
boundary  conditions  on  the  object  and  numerically  solving  the 
integral  equations  that  describe  the  surface  current.  Recall 
the  magnetic  field  integral  equation  given  by  (1). 
Simulations  produced  by  Morgan's  time-domain  thin  wire 
integral  equation  computer  program  [12]  were  used  to  evaluate 
the  pole  extraction  algorithm.  The  excitation  waveform  used 
is  the  double  Gaussian  pulse  depicted  in  Figure  9.  This  pulse 
is  a  wide  Gaussian  pulse  with  a  ten  percent  width  of  0.3 
nanoseconds  subtracted  from  a  narrow  Gaussian  pulse  with  a  ten 
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nanoseconds  subtracted  from  a  narrow  Gaussian  pulse  with  a  ten 
percent  width  of  0.15  nanoseconds. 

Figures  10  through  13  depict  back  scattering 
response  of  a  0.1  meter  length  thin-wire,  having  a  radius  of 
0.00118  meter,  computed  at  various  incident  aspects,  ranging 
from  thirty  degrees  to  ninety  degrees.  The  laboratory 
arrangement  for  actual  measurements  simulated  by  Morgan's 
program  is  described  in  [17]  .  Ninety  degrees  represents  a 
broadside  aspect,  while  thirty  degrees  represents  the  incident 
plane  wave  having  nearly  grazing  incidence  on  the  wire.  The 
poles  extracted  at  each  of  the  four  aspect  angles  are  plotted 
in  Figure  14.  In  this  figure,  and  those  that  follow  which 
depict  extracted  poles,  the  signal  poles  lie  in  or  on  the  unit 
circle,  and  the  noise  poles  lie  outside. 

The  results  obtained  with  this  rigorous  numerical 
computation  demonstrate  the  aspect  independence  of  the 
extracted  poles  using  the  Kumaresan-Tuf ts  method.  Note  that 
only  half  of  the  poles  were  obtained  for  broadside 
illumination;  two  even-numbered  poles  can  easily  be  seen 
outside  the  unit  circle.  This  results  because  of  the  physical 
symmetry  of  both  the  wire  and  the  incident  field,  thus 
precluding  excitation  of  odd-symmetric  modal  currents  and 
their  associated  natural  resonances. 

Figure  15  exemplifies  the  computed  back-scattering 
response  of  the  0.1  meter  thin  wire  corrupted  artificially 
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Figure  10.   Integral  Equation  Thin  Wire  Scattering,  30  Degree 
Aspect 


32 


X103  45  degree  backscattering 

4 


3  - 


2  - 


1  - 


-1  - 


-2- 


-3  - 


-5 


0  0.5  1  1.5  2  2.5  3  3.5  4  4.5  5 

Time   (nanoseconds) 


Figure  11.   Integral  Equation  Thin  Wire  Scattering,  45  Degree 
Aspect 
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Figure  12.   Integral  Eguation  Thin 
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Figure  14.   Kumaresan-Tufts  Poles,  Noiseless  Thin  Wire  Data 
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Figure  15.   Integral  Eguation  Thin  Wire  Scattering,  20.0  dB 
SNR,  4  5  Degree  Aspect 
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with  noise  at  a  20.0  dB  SNR.  Figure  16  shows  the  poles 
extracted  at  each  of  the  four  angles  of  incidence  used 
previously  in  Figure  14.  Poles  of  Figure  14  at  90°  are  now 
missing  in  Figure  16,  and  only  the  first  three  low  frequency 
poles  are  tightly  grouped.  The  loss  of  high  frequency  poles 
is  expected  because  these  have  the  highest  damping  and  thus 
lose  their  energy  at  the  fastest  rate.  Further  comparison 
between  results  computed  at  20.0  dB  SNR  and  infinite  SNR  are 
offered,  angle  by  angle,  in  Figures  17  through  20. 

One  additional  test  of  the  computed  thin  wire 
scattering  was  conducted  at  a  7.0  dB  SNR.  The  corrupted 
waveforms  are  exemplified  by  Figure  21;  the  extracted  poles 
are  shown  in  Figure  22.  The  number  of  poles  obtained  has 
decreased  with  respect  to  the  number  obtained  at  20.0  dB  SNR. 
The  grouping  of  the  clusters  has  also  expanded.  Angle  by 
angle  comparisons  are  again  offered  in  Figures  23  through  26. 
c.   Scale  Model  Measurements 

The  transient  scattering  measurements  of  scale 
models  used  for  evaluation  in  this  section  were  made  by  Walsh 
using  the  anechoic  chamber  of  the  Transient  Electromagnetic 
Scattering  Laboratory  at  the  Naval  Postgraduate  School.  The 
entire  measurement  process  and  laboratory  setup  are  described 
in  detail  in  [17] . 


38 


Extracted  poles 


c 
'5b 


Real  z 


Figure  16.   Kumaresan-Tufts  Poles,  20.0  dB  SNR 
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Figure  17.   Integral  Eguation  Thin  Wire  Comparison,  Noiseless 
vs.  20.0  dB  SNR,  30  Degree  Aspect 
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Figure  18.   Integral  Equation  Thin  Wire  Comparison,  Noiseless 
vs.  2  0.0  dB  SNR,  4  5  Degree  Aspect 
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Figure  19.   Integral  Eguation  Thin  Wire  Comparison,  Noiseless 


vs.  20.0  dB  SNR,  60  Degree  Aspect 
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Figure  20.   Integral  Equation  Thin  Wire  Comparison,  Noiseless 
vs.  2  0.0  dB  SNR,  9  0  Degree  Aspect 
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Figure  21.   Integral  Eguation  Thin  Wire  Scattering,  7.0  dB 
SNR,  4  5  Degree  Aspect 
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Figure  22.   Kumaresan-Tufts  Poles,  7.0  dB  SNR 
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Figure  23.   Integral  Equation  Thin  Wire  Comparison,  Noiseless 
vs.  7.0  dB  SNR,  30  Degree  Aspect 
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Figure  24.   Integral  Equation  Thin  Wire  Comparison,  Noiseless 
vs.  7.0  dB  SNR,  45  Degree  Aspect 
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Figure  25.   Integral  Equation  Thin  Wire  Comparison,  Noiseless 
vs.  7.0  dB  SNR,  60  Degree  Aspect 
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Figure  26.   Integral  Equation  Thin  Wire  Comparison,  Noiseless 
vs.  7.0  dB  SNR,  90  Degree  Aspect 
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1.   Wire  Targets 

The  thin  wire  measurements  were  obtained  from 
the  scattering  response  of  a  0.1  meter  length  thin  wire  having 
radius  0.00118  meter.  Recall  that  these  are  the  same 
dimensions  as  the  wire  whose  computed  response  was  processed 
in  the  previous  section.  The  measurements  at  each  of  four 
incident  aspects  are  shown  in  Figures  27  through  30. 

The  poles  extracted  from  the  four  measurements 
are  depicted  in  Figure  31.  As  before  in  the  computed  noisy 
data,  tight  clusters  occur  only  at  the  lowest  frequencies. 
The  poles  in  these  tight  clusters  are  those  which  are 
measurably  present  at  various  aspects.  The  poles  extracted 
at  higher  frequencies  are  those  which  possessed  sufficient 
measurable  energy  at  the  given  aspect.  Figure  32  depicts  the 
comparison  between  poles  extracted  from  the  measured  and 
computed  signals.  Again,  the  closest  agreement  between  the 
two  sets  of  poles  occurs  at  the  lowest  frequences. 
2.   Aircraft  Models 

Plastic  1/72  scale  aircraft  models,  coated  with 
silver,  were  used  for  transient  scattering  measurements. 
Representative  scattering  signatures  of  two  aircraft  targets, 
measured  at  six  different  aspects,  are  shown  in  Figures  33 
through  36. 

The  results  of  pole  extraction  in  target  1  are 
shown  for  a  total  of  six  different  aspects  in  Figures  37  and 
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Figure  27.   Measured  Thin  Wire  Scattering,  30  Degree  Aspect 
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Figure  28.   Measured  Thin  Wire  Scattering,  45  Degree  Aspect 
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Figure  29.   Measured  Thin  Wire  Scattering,  60  Degree  Aspect 
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Figure  31.   Kumaresan-Tufts  Poles,  Measured  Thin 


Wire 
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Figure   32.  Thin   Wire   Comparison,   Measured   vs.   Integral 
Equation 
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Figure  33.   Target  1  Scattering,  30  Degrees  from  Nose  on 
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Figure  34.   Target  1  Scattering,  Nose  on 
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Figure  35.   Target  2  Scattering,  30  Degrees  from  Nose  on 
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Figure  36.   Target  2  Scattering,  Nose 
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Figure  37.   Kumaresan-Tufts  Poles,  Target  1 
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38.  The  poles  extracted  at  all  six  aspects  are  shown  in 
Figure  39.  Only  one  clearly  discernible  cluster  is  present 
in  each  of  the  three  figures.  At  higher  frequencies,  no 
useful  information  is  imparted  by  the  data.  Results  of 
similar,  though  slightly  improved  quality,  were  obtained  from 
target  2.  These  results  are  presented  in  Figures  40  through 
42  in  the  format  of  Figures  37  through  39  respectively. 

Although  the  Kumaresan-Tuf ts  algorithm  is 
capable  of  extracting  low  frequency  poles  acceptably,  the 
inconsistent  results  at  higher  frequences  reveals  the  inherent 
weakness  in  an  algorithm  capable  of  processing  only  the  late- 
time  portion  of  a  target's  radar  response. 

A  side-by-side  comparison  of  poles  obtained  from 
both  aircraft  by  both  the  Kumaresan-Tuf  ts  method  and  the 
Cadzow-Solomon  method  is  presented  at  the  end  of  the  chapter 
to  illustrate  the  gains  afforded  by  processing  the  early— time. 

C.   CADZOW-SOLOMON  ALGORITHM 

Recall  from  the  results  depicted  in  Figure  8  that  a  late 
transition  to  late-time,  and  the  consequent  reduction  of 
signal  power,  caused  complete  breakdown  of  the  Kumerasan-Tuf ts 
algorithm.  The  Cadzow-Solomon  algorithm  addresses  this 
shortcoming  by  processing  the  signal  at  the  instantaneous 
onset  of  early-time.  Thus,  the  Cadzow-Solomon  algorithm  is 
capable  of  processing  the  earliest  response  of  a  target  to 
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Figure  38.   Kumaresan-Tufts  Poles,  Target  1 
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Figure  39.   Kumaresan-Tufts  Poles,  Target  1 
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Figure  40.   Kumaresan-Tufts  Poles,  Target  2 
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Figure  41.   Kumaresan-Tufts  Poles,  Target  2 
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Figure  42.   Kumaresan-Tufts  Poles,  Target  2 
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electromagnetic  excitation  where  the  response  has  the  greatest 
magnitude . 

1.  Applicability 

The  early-time  portion  of  a  target's  scattered  field 
occurs  as  long  as  there  is  a  driven  portion  of  the  total 
field.  Once  the  field  no  longer  contains  a  scattered  response 
due,  in  part,  to  the  incident  excitation  at  points  on  the 
object,  early-time  ceases  and  late— time  begins.  Hence,  the 
Cadzow-Solomon  models  both  the  system's  input  and  output,  and 
equivalently ,  the  poles  and  zeros  of  the  system  transfer 
function. 

2.  Equations 

The  Cadzow-Solomon  algorithm  extends  the  auto- 
regressive  equation  (9),  used  in  Prony's  method,  to  the  more 
general  autoregressive  moving  average  (ARMA)  equation 


y  =Sbiy„.,+  Sa1xn 


1  =  1 


n-1 


(23) 


1=0 


where  the  second  summation  term  models  the  excitation  to  the 
system. 

A  set  of  M  such  equations  in  matrix  form  is  given  by 


YKD-1   X0  • 


K-r-Y*'+H-2x"->- '**»+"-> . 


L  a* 


|_  aKD+M-'  _1 


0   _I 


(24) 
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As  in  the  Kumaresan-Tuf ts  method,  M  is  selected  to  be  greater 
than  the  column  dimension  of  the  data  matrix  which  is 
KD+KN+1  . 

3.   Excess  Poles  and  Noise  Removal 

The  Cadzow-Solomon  method  used  in  this  thesis  is  a 
modification  which  incorporates  the  non-causal  arrangement  of 
the  system  equations  used  by  Kumaresan-Tuf ts .  This 
modification  was  first  discussed  by  Norton  in  [16]  .  The 
Kumaresan  approach  of  overestimating  the  system  order  can  be 
used  as  before  in  a  non-causal  model  to  constrain  the  noise 
poles  inside  the  unit  circle,  while  SVD  forces  the  signal 
poles  outside  the  unit  circle. 

Since  the  input  waveform  is  known,  its  order  can  be 
almost  exactly  determined.  In  all  the  work  of  this  thesis, 
the  input  waveform  used  is  the  double  Gaussian  depicted  in 
Figure  14.  Approximately  25  samples  defining  this  pulse  of 
0.5  nanoseconds  duration  makes  KN  equal  25  in  equation  (23). 
Since  the  input  is  causal,  the  signal  zeros  fall  inside  the 
unit  circle  where  they  cannot  be  easily  segregated  from 
similarly  located  noise  poles.  However,  the  signal  zeros 
impart  no  information  about  the  target  and  need  not  be 
extracted.  The  inclusion  of  the  input  in  the  data  matrix  is 
nevertheless  vital  to  the  model  of  the  system  and  the  accurate 
determination  of  the  signal  poles. 
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The  ARMA  equation  of  (23)  can  be  modified  to  obtain 


y  =  Y  b'.y     +7  a,xn 

n  1=1     KD  n  '  l   1=0 


(25) 


The  recursive  portion  of  (25)  is  now  in  a  non-causal  form 
similar  to  expression  (12)  .  A  set  of  M  such  equations  in 
matrix  form  is  given  by 


KN  +  1 


KN  +  KD 


Ar.   ■  §  •    «v 


y   •  •  •  y     xm- 


Or,  in  matrix  notation 


.  .  x 


KN+M-1 


_     *:o     J 


.  ^V""' 


(26) 


[Dy*][-5-]=y    «> 


ere 


[Dyx]  =  [Dy:Dx] 


(27) 


4.  Singular  Value  Decomposition 

Like  the  system  equations  of  the  Kumaresan-Tuf ts 
model,  the  system  equations  in  (26)  are  processed  using 
singular  value  decomposition.  The  coefficient  vector  is  again 
the  minimum-norm  solution,  which  constrains  the  extraneous 
poles  and  extraneous  zeros  to  be  inside  the  unit  circle. 

5.  Bias  Compensation  in  the  Cadzow-Solomon  Formulation 
By  compensating  the  eigenvalues  of  the  i  matrix  in 

(16),  the  performance  of  the  Kumaresan-Tuf ts  algorithm  is 
significantly  improved  in  the  presence  of  noise.  Cadzow- 
Solomon  have  shown  [5]  that  if  the  actual  orders  K^  and  KN 
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are  overestimated  to  be  kd  and  kn  ,  min  (kd~kd  >  kn~kn)  singular 
values  are  zero  in  noiseless  data.  Since  the  input  data  is 
known,  the  eigenvalues  of  the  data  matrix  may  be  compensated 
in  the  same  manner  as  in  the  Kumaresan-Tuf ts  algorithm  for 
noiseless  data. 

To  understand  the  compensation  required  in  noisy  data, 
an  analysis  of  additive  noise  is  required.  As  given  by  Norton 
[16] ,  if  the  input  data  noise  is  W(  and  the  output  data  noise 
is  v,  r    the  data  matrix  may  be  modeled  as 


where 


and 


[Dyx]  =  [Dy:Dx]=Syx+Ny> 


[Nyx]  =  [Ny:Nx] 


(28) 


(29) 


Nx= 


WM  •  *  *  WM  + 


M  +  Kr 


Ny= 


.V 


M  +  K 


D  _l 


(30) 


The  expected  value  of  D  DT     is  then 

F  yx  yx 


E[DvxDTyx]=SyxST+E[NyxNTyx] 


(31) 


If  the  input  and  output  noise  variances  are  not  equal,  the 
eigenvalue  shifting  theorem  used  in  Kumaresan-Tuf ts  cannot  be 
used  to  analytically  predict  the  requisite  eigenvalue 
compensation  of  D„YDlr.-    Nevertheless,  when  the  input  and 

yx  yx 

output  variances  were  assumed  equal,  and  eigenvalue 
compensation  similar  to  that  used  in  Kumaresan-Tuf ts  was 
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performed,  the  results  were  consistently  superior  to  those 
obtained  without  compensation.  Therefore,  the  results  of 
Cadzow-Solomon  signal  processing  presented  in  this  thesis  were 
obtained  using  eigenvalue  compensation  and  the  assumption  of 
equal  noise  variance. 
6.   Performance 

The  Cadzow-Solomon  algorithm  was  programmed  in  Fortran 
and  tested  on  the  same  data  used  for  evaluating  the  Kumaresan- 
Tufts  algorithm.  Note  that  the  Cadzow-Solomon  algorithm 
can  use  the  early-time  portion  of  the  data  that  the  Kumaresan- 
Tufts  algorithm  can  not  use.  The  program  appears  in 
Appendix  B. 

a.   Synthetically  Generated  Data 

The  starting  point  for  evaluating  the  performance 
of   the   Cadzow-Solomon   algorithm   was   with   synthetically 
generated  data  of  the  form  given  by  (8)  plus  the  addition  of 
input  data  required  to  model  early  time  data. 
1.   Noise  Performance 

The  algorithm  was  evaluated  at  various  signal- 
to-noise  ratios,  ranging  from  90.0  dB  to  7.0  dB.  Figure  43 
shows  the  signal  produced  by  two  s-plane  poles  at  90.0  dB , 
with  a  late-time  beginning  at  10.0  nanoseconds.  Figures  44 
through  48  depict  the  poles  extracted  from  this  signal  at  the 
different  signal-to-noise  ratios. 
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Figure  43.   Signal  Containing  Two  S-Plane  Poles, 
90.0  dB  SNR 
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Figure  44.   Cadzow-Solomon  Poles,  Synthetic  Data,  90.0  dB  SNR 
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Figure  45.   Cadzow-Solomon  Poles,  Synthetic  Data,  30.0  dB  SNR 
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Figure  46.   Cadzow-Solomon  Poles,  Synthetic  Data,  20.0  dB  SNR 
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Figure  47.   Cadzow-Solomon  Poles,  Synthetic  Data,  10.0  dB  SNR 
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Figure  48.   Cadzow-Solomon  Poles,  Synthetic  Data,  7.0  dB  SNR 
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The  figures  chart  the  steady  degradation  of 
the  algorithm's  performance  with  the  increase  of  noise.  At 
30.0  dB ,  the  location  of  the  low  frequency  pole  is  already 
slightly  displaced.  More  significant  is  the  location  of  one 
of  the  extracted  poles  in  the  noise  signal  space.  At  20.0  dB , 
the  low  frequency  pole  is  located  in  some  trials  on  the  real 
axis.  At  10.0  dB ,  all  the  extractions  are  located  on  the  real 
axis  and  at  7.0  dB  their  locations  there  are  dispersed.  The 
extraction  of  the  higher  frequency  pole  is 
uncharacteristically  more  accurate  than  that  of  the  low 
frequency  pole.  Even  at  7.0  dB ,  the  high  frequency  pole  is 
located  with  excellent  accuracy.  The  location  of  the  low 
frequency  pole  near  the  real  axis  was  chosen  deliberately  to 
illustrate  the  difficulty  in  resolving  the  slight  frequency 
difference  between  the  true  pole  and  a  noise  pole  located  on 
the  real  axis.  Also,  fewer  points  were  processed  using  the 
Cadzow-Solomon  method  than  were  processed  using  the  Kumaresan- 
Tufts  method,  since  the  largest  data  matrix  allowed  by  the 
programs  in  Appendices  A  and  B  contain  fewer  data  points  in 
the  Cadzow-Solomon  data  matrix  than  in  the  Kumaresan-Tuf ts 
data  matrix.  The  results  demonstrate  the  need  to  process  a 
substantial  number  of  points  in  order  to  accurately  extract 
low  frequency  poles. 
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b.  Thin  Wire  Integral  Equation  Generated  Data 

The  performance  of  the  Cadzow-Solomon  algorithm 
was  evaluated  using  the  same  set  of  data  tested  by  the 
Kumaresan-Tuf ts  algorithm.  The  results  are  presented  in 
Figure  49.  Tight  clusters  appear  at  frequencies  higher  than 
those  obtained  with  the  Kumaresan-Tuf ts  algorithm.  Figure  50 
depicts  the  poles  extracted  from  the  same  signal  at  a  20.0  dB 
SNR.  The  clustering  at  this  SNR  is  comparable  to  the  results 
obtained  by  the  Kumaresan-Tuf ts  method  with  the  noiseless 
data.  Further  angle-by-angle  comparisons  of  the  poles 
extracted  from  the  noiseless  data  and  the  20.0  dB  data  are 
depicted  in  Figures  51  through  54.  Note  the  small  number  of 
poles  in  Figure  54  due  to  the  unexcited  odd-symmetric  poles 
at  90°  aspect. 

One  further  test  was  conducted  on  computed  data 
at  a  7.0  dB  SNR.  The  results  are  depicted  in  Figure  55.  Even 
at  7.0  dB ,  discernible  clusters  are  present.  Angle-by-angle 
comparisons  of  the  poles  obtained  in  7.0  dB  data  and  those 
obtained  in  noiseless  data  are  presented  in  Figures  56  through 
59. 

c.  Scale  Models 

The  same  scale  models  used  to  evaluate  the 
Kumaresan-Tuf ts  algorithm  were  used  to  evaluate  the  Cadzow- 
Solomon  algorithm. 


80 


Extracted  poles 


c 
'5b 

eo 

E 


0.5  - 


0^ 


-0.5  - 


-1  - 


1                  i                   '       "  7 1 1 

& 

+ 
+         X 

+     / 
+     / 

+ 

-^                               o 

^*                                                             X 

* 

*  \ 

+ 
+ 

K 

+ 

9 

+ 

+  30  degrees 

\  ° 

\ 

m\ 

— 

\ 

x  45  degrees 
*  60  degrees 
o  90  degrees 

*/ 

/  ° 

•  / 

+     V 

*      \ 

*      \ 

+         ^v 

+ 

xf          ^^*^" 

*  / 

/  0 

: 

-0.5 


0.5 


Real  z 


Figure  49.   Cadzow-Solomon  Poles,  Noiseless  Thin  Wire  Data 
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Figure  50.   Cadzow-Solomon  Poles,  20.0  dB  SNR 
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Figure  51.   Integral  Equation  Thin  Wire  Comparison,  Noiseless 
vs.  20.0  dB  SNR,  30  Degree  Aspect 
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Figure  52.   Integral  Equation  Thin  Wire  Comparison,  Noiseless 
vs.  2  0.0  dB  SNR,  4  5  Degree  Aspect 
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Figure  53.   Integral  Equation  Thin  Wire  Comparison,  Noiseless, 
vs.  20.0  dB  SNR,  60  Degree  Aspect 
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Figure  54.   Integral  Equation  Thin  Wire  Comparison,  Noiseless 
vs.  20.0  dB  SNR,  90  Degree  Aspect 
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Figure  55.   Cadzow-Solomon  Poles,  7.0  dB  SNR 
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Figure  56.   Integral  Equation  Thin  Wire  Comparison,  Noiseless 
vs.  7.0  dB  SNR,  3  0  Degree  Aspect 


88 


Extracted  poles 


1  - 


0.5  - 


C3 
C 

'ob 
rt 

E 


0- 


-0.5 


-1  - 


1                                          1                                          1 

i 

o             -—— — 

°     °  >^ 

+      °      /< 

:  +         / 

+ 

o 

+  ^v 

X 

+    :          / 

O:  / 

0       j 

A 

i 

■h-e 

0 

0  7  dB  SNR 

0    \ 

+          \. 

i 
V 

9  / 

i  +       \ 

+       O        N. 

+                 \. 

0         \. 

O                                           \^^                                                   + 

+              ^ — ^ 

+ 

V 

O 

!        i     ■       • 
i               i               i 

-0.5 


0.5 


Real  z 


Figure  57.   Integral  Equation  Thin  Wire  Comparison,  Noiseless 
vs.  7.0  dB  SNR,  4  5  Degree  Aspect 
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Fiaure  58    Integral  Equation  Thin  Wirecompa 
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Fiaure  59    Integral  Equation  Thin  Wire  Comparison,  Noiseless 
y       '   vs.  7.0  dB  SNR,  90  Degree  Aspect 
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1.  Wire  Targets 

Figure  60  depicts  the  poles  extracted  from 
measurements  of  a  0.1  meter  wire.  Three  tight  clusters  appear 
at  the  lowest  frequencies  and  at  the  highest  frequencies.  The 
poles  in  between  can  not  be  easily  discriminated.  The 
dispersion  of  these  poles  is  apparently  due  to  the  aspect 
dependence  of  their  measurable  power.  In  other  words,  these 
poles  are  excited  more  at  some  aspects  then  at  others. 

Figure  61  depicts  the  comparison  between  poles 
extracted  from  computed  data  and  measured  data.  As  in  Figure 
60,  close  agreement  exists  at  the  highest  and  lowest 
frequencies.  The  results  are  much  more  favorable  than  those 
similarly  obtained  by  the  Kumaresan-Tuf ts  algorithm. 

2.  Model  Aircraft 

Figures  62  through  64  depict  poles  extracted 
from  aircraft  target  1.  As  in  the  Kumaresan-Tuf ts  testing, 
the  Cadzow-Solomon  testing  was  conducted  at  six  different 
aspects.  Results  for  target  2  are  depicted  in  Figures  65 
through  67.  The  results  of  both  targets  show  clearly  defined 
clusters.  The  first  two  clusters  of  target  2  are 
exceptionally  tight.  However,  the  mid-frequency  clusters  of 
target  2  are  not  as  clearly  formed  as  those  of  target  1. 

Comparisons  of  poles  obtained  with  each  method 
for  target  1  and  2  are  depicted  in  Figure  68  and  69 
respectively.   These  two  figures  graphically  depict  the  clear 
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Figure  60.   Cadzow-Solomon  Poles,  Measured  Thin  Wire 
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Figure   61.  Thin  Wire   Comparison,   Measured  vs.   Integral 
Equation 
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Figure  62.   Cadzow-Solomon  Poles  Target  1,  Three  Aspects 
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Figure  63.   Cadzow-Solomon  Poles  Target  1,  Three  Aspects 
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Figure  64.   Cadzow-Solomon  Poles  Target  1,  All  Six  Aspects 
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Figure  65.   Cadzow- Solomon  Poles  Target  2,  Three  Aspects 
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Figure  66.   Cadzow-Solomon  Poles  Target  2,  Three  Aspects 
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Figure  67.   Cadzow-Solomon  Poles  Target  2,  All  Six  Targets 
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Figure  68.   Pole  Comparisons,  Target  1,  All  Six  Targets 
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Figure  69.   Pole  Comparisons,  Target  2,  All  Six  Aspects 
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superiority   of   the   Cadzow-Solomon   algorithm   over   the 
Kumaresan-Tuf ts  algorithm. 

In  order  to  obtain  an  initial  indication  of 
the  possibility  for  target  classification  through  pole 
extraction,  nose-on  measurements  of  two  additional  aircraft 
models  were  made,  processed  and  compared  with  the  results  of 
targets  1  and  2.  The  nose-on  measurements  of  targets  3  and 
4  appear  in  Figures  70  and  71  respectively.  A  comparison  plot 
of  poles  extracted  from  each  of  the  four  targets  is  depicted 
in  Figure  72.  Each  of  the  four  aircraft  measured  are  fighters 
of  similar  size  and  shape  (see  Table  1).  The  poles  for  each 
target  are  sufficiently  different  in  this  single  measurement 
to  identify  each  aircraft  individually.  However,  some  of  the 
poles  are  arranged  in  clusters  which  appear  with  a  harmonic 
pattern  similar  to  that  obtained  for  either  of  the  first  two 
aircraft  at  various  aspects.  In  order  to  more  fully  assess 
the  target  classification  capability  of  pole  extraction, 
several  measurements  should  be  made  of  a  given  aircraft  model. 
A  plot  of  the  poles  extracted  from  each  of  these  measurements 
would  form  clusters  at  the  locations  of  the  true  poles.  The 
centroid  of  each  of  these  clusters  would  then  be  compared 
against  the  centroid  poles  similarly  obtained  from  other 
aircraft.  Although  several  poles  of  different  aircraft  might 
be  similar,  the  set  of  poles  belonging  to  an  aircraft  could 
form  the  basis  for  classification  if  that  set  was  unique  among 
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Figure  72.   Cadzow-  Solomon  Pole  Comparisons,  4  Targets,  Nose-on 
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the  sets  belonging  to  all  other  measured  aircraft.  The  results 
in  Figure  72  demonstrate  the  possibility  of  using  the  Cadzow- 
Solomon  pole  extraction  algorithm  to  aid  in  the  classification 
of  aircraft,  perhaps  by  use  of  the  extracted  poles  in 
constructing  annihilation  filters. 

TABLE  1.   FULL  SIZE  DIMENSIONS  OF  TARGETS  RECORDED 


Target  number 

L 

2 

3 

i 

4 

Overall  length 

12 

.20 

15 

.03 

16 

.94 

16 

.00 

(meters ) 

Overall  height 

3 

.35 

5 

.09 

4 

.51 

4 

.80 

(meters ) 

Wingspan 

10 

.96 

10 

.00 

11 

.43 

13 

.95 

(meters ) 

Tailplane  span     Unknown      5.58        6.92         5.75 
(meters ) 
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III.   SUMMARIES  AND  CONCLUSIONS 

In  this  chapter,  a  step-by-step  guide  through  each 
algorithm  is  presented.  At  each  step,  techniques  and  lessons 
learned  are  discussed  together  with  general  observations. 
Conclusions  are  presented  at  the  end  of  the  chapter. 

A.   KUMARESAN-TUFTS 

The  first  step  in  processing  a  signal  with  the  Kumaresan- 
Tufts  algorithm  is  to  determine  the  beginning  of  early-time. 
The  objective  is  to  pick  the  earliest  possible  starting  point 
without  entering  into  the  latter  part  of  early-time.  If  the 
starting  point  for  processing  is  improperly  chosen  to  include 
the  early-time,  the  results  will  be  completely  unreliable 
since  the  signal  no  longer  satisfies  the  late  time  model.  If 
the  starting  point  is  chosen  too  late,  the  signal  may  not  be 
sufficiently  strong  in  the  presence  of  measurement  noise. 
Since  the  signal  is  the  sum  of  exponentially  damped  sinusoids, 
the  optimum  starting  point  is  at  the  precise  instant  of 
transaction  into  late-time.  The  key  to  determining  the 
beginning  of  late-time  is  in  determining  the  beginning  of 
early  time.  Determining  the  first  response  of  the  target  to 
excitation  cannot  usually  be  done  by  a  simple  visual 
inspection  of  measurement  data.   Unless  the  exact  distance  to 
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the  target  is  known,  the  most  accurate  method  attempted  by  the 
author  for  determining  the  beginning  of  early— time  is  to 
process  the  signal  using  the  Cadzow-Solomon  algorithm.  This 
is  discussed  in  the  next  section.  However,  the  reliance  of 
the  Kumaresan-Tuf ts  algorithm  on  information  provided  by  the 
Cadzow-Solomon  algorithm  is  an  obvious  disadvantage  of  the 
former  method. 

Once  the  starting  point  for  processing  has  been  selected, 
the  next  step  is  to  determine  the  dimensions  of  the  data 
matrix  and,  consequently,  the  number  of  points  in  the  signal 
to  be  processed.  In  trials  conducted  on  noiseless  synthetic 
data,  the  accuracy  of  pole  extraction  increased  steadily  with 
the  increase  in  the  data  matrix  dimensions.  These  trials  were 
conducted  up  to  the  limit  of  the  array  dimensions  defined  in 
the  computer  program  of  Appendix  A.  The  number  of  points 
processed  in  measurement  data  should  be  as  large  as  possible, 
while  still  meeting  the  following  two  constraints.  First, 
incorporate  as  many  cycles  of  the  data  as  possible.  Usually, 
visual  inspection  of  the  data  reveals  a  repeating  pattern 
which  should  be  entirely  incorporated  into  the  window  of 
points  to  be  processed.  When  only  portions  of  these  patterns 
are  selected,  a  disproportionate  weighting  tends  to  be  placed 
on  certain  poles.  Second,  signal  portions  late  in  the 
response  which  are  no  longer  distinguishable  in  the  presence 
of  measurement  noise  should  not  be  selected. 
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The  final  step  involves  determining  the  number  of  true 
poles  in  the  system.  The  following  approach  has  proven  to  be 
the  most  successful.  First,  process  the  signal  without  any 
eigenvalue  compensation  to  establish  an  upper  bound  on  the 
order  of  the  system.  In  most  cases,  the  number  of  poles 
outside  the  unit  circle  will  be  less  than  the  overestimated 
order  of  the  system.  If  not,  increase  the  row  dimension  of 
the  data  matrix  in  order  to  increase  the  estimated  order  of 
the  system,  and  repeat.  When  the  number  of  poles  is  less  than 
the  estimated  order  of  the  system,  then  one  should  gradually 
increase  the  number  of  eigenvalues  compensated  in  successive 
trials,  while  closely  observing  the  effects  induced  on  the 
poles  outside  the  unit  circle.  As  the  number  of  eigenvalues 
compensated  is  steadily  increased,  noise  poles  and  weak  signal 
poles  will  move  inside  the  unit  circle.  The  programs  in 
Appendix  A  and  B  allow  the  user  to  compare  the  results  of 
successive  trials,  by  generating  overlays  for  each  plot.  If 
N  poles  are  in  the  signal  space,  at  least  the  first  N 
eigenvalues  must  not  be  compensated,  or  true  poles  may  be 
lost.  As  the  actual  order  of  the  system  is  approached  by 
compensation,  the  user  will  notice  an  orderly,  even 
arrangement  assumed  by  the  noise  poles.  If  certain  poles 
still  remain  suspect  after  compensation,  vary  slightly  the 
other   parameters,   such   as   the   starting   point   and   the 
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dimensions  of  the  data  matrix.  Generally,  only  true  signal 
poles  will  repeatedly  assert  themselves  under  varying 
parameters . 

B.   CADZOW-SOLOMON 

The  techniques  and  general  observations  offered  in  the 
preceding  section  apply  equally  to  the  Cadzow-Solomon 
algorithm.  An  important  consideration  in  this  method,  not 
discussed  above,  is  the  selection  of  the  beginning  of  early- 
time.  Candidates  for  a  starting  point  are  usually  at  or  near 
zero  crossings  within  approximately  thirty  points  of  the 
object's  first  definite  response  to  electromagnetic 
excitation.  Begin  processing  at  the  chosen  point  while 
varying  parameters  in  successive  trials.  Select  the  point 
whose  successive  results  are  the  most  consistent  under  varying 
parameters . 

The  selection  of  the  starting  point  for  beginning  of 
early-time  can  be  very  critical.  For  example,  not  a  single 
pole  could  be  extracted  in  one  trial  wherein  the  starting 
point  occurred  only  ten  points  after  the  actual  starting 
point.  Additionally,  in  most  cases  observed,  the  late-time 
start  given  by  the  selected  early-time  occurred  within  less 
than  two  points  from  a  zero  crossing.  If  this  observation 
proves  to  be  generally  true  in  later  research,  it  may  serve 
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as  a  way  to  check  the  starting  point  selected  for  one 
algorithm  in  terms  of  the  other. 

C.   CONCLUSIONS 

Both  the  Kumaresan-Tuf ts  and  the  Cadzow-Solomon  algorithms 
can  effectively  extract  poles  from  the  scattering  response  of 
a  radar  target.  Because  both  algorithms  obtain  a  least- 
squares  solution  to  the  system  model,  both  perform  acceptably 
in  the  presence  of  noise.  Although  eigenvalue  compensation 
is  not  analytically  justified  in  the  Cadzow-Solomon  algorithm, 
the  results  obtained  through  eigenvalue  compensation  in  this 
method  were  generally  superior  to  those  similarly  obtained  in 
the  Kumaresan-Tuf ts  method.  The  results  demonstrated  the 
inherent  advantages  of  an  algorithm  capable  of  processing  a 
target's  strongest  response  in  the  early  time.  The  Kumaresan- 
Tufts  method  compared  favorably  with  the  Cadzow-Solomon  only 
in  responses  with  a  long  late— time. 
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APPENDIX  A.   THE  KUMARE SAN-TUFTS  POLE  EXTRACTION  ALGORTITHM 


The  following  program  implements  the  Kumaresan-Tuf ts 
algorithm  as  described  in  Chapter  2  of  this  thesis.  The 
program  is  written  in  Fortran  77.  The  SVD  and  root-finding 
subroutines  called  by  this  program  are  found  in  the  EISPACK 
library  [18] .  The  SVD  subroutine  is  a  translation  from  ALGOL 
as  given  in  [19] .  The  matrix  multiplication  and  graphics 
subroutines,  also  called  by  this  program,  are  found  in 
Appendix  C  and  D  respectively. 


INTEGER  IERR,Kd,M,MN,MAGPCL,NSTRTPT,DELTAY 

INTEGER  IER,NCAUS,NMENU,L/1/ 

Dm3GER*2  KdPLT 

REAL*8  A(70,70) ,W(70) ,U(70,70) ,V(70,70) ,RV1(70) 

REAL*8  VS (70,70) ,UT(70,70) ,AINV(70,70) ,X(70) 

REAL*8  XP(70) ,B(70) ,SIGMA(70,70) ,SIG(70,70) 

REAL*8  COF(70),RCCTR(70),ROOTI(70) 

REAL*8  D(1024) ,AVG,MACHEP/1.0E-16/,Dy(140) 

C0MPLEX*16  S(70) 

LOGICAL  MATU/. TRUE. /,MATV/.TRUE./, CAUSAL/. TRUE./, LONG/. TRUE./ 

LOGICAL  DSET/. FALSE. /,NUFILE/. TRUE./ 

CHARACTER  TITLE*16,HEAI^*64,YI^l,DC*l,TITTJRn6,TITLEI*16 

CHARACTER  TrTL*16 


C    Enter  parameters  for  processing 

11         IF  (DSET)  CLOSE  (10) 
N0VERLAY=0 
CFEN(10,FILE=,PLC'r) 


IF  (DSET 

WRITE  (* 

WRITE  (* 

WRITE  (* 

WRITE  (* 

WRITE  (* 

WRITE  (* 

WRITE  (* 

WRITE  (* 


GO  TO  85 


Welcome  to  signal  processing  using  the' 
Kumaresan-Tuf ts  method' 

Do  you  want  ' 

1.  The  long  version  for  beginners' 

2.  The  short  version  for  pros' 
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15   WRITE  (*,*)  'Please  enter  1  or  2  ' 
READ  (*,*)  N 
IF  (N  .EQ.  1)  THEN 
LO¥G=.TRUE. 

ELSEIF  (N  .EQ.  2)  THEN 
I£NG=. FALSE. 
ELSE 
GOTO  15 
ENDIF 

WRITE  (*,*)  'Session  will  begin  with  entry  of  parameters  needed  fo+r  processing' 
WRITE  (*,*) 


16 


10 


WRITE  (*,*) 

Do  you  want  to  enter  parameters 

from' 

WRITE  (*,*) 

i  i 

WRITE  (*,*) 

1.  The  keyboard' 

WRITE  (*,*) 

2.  A  previously  created  file  of 

parameters 

WRITE  (*,*) 

i 

WRITE  (*,*) 

Please  enter  1  or  2  ' 

READ  (*,*)  N 

IF  (N  .EQ.  1) 

THEN 

GO  TO  1 

ELSEIF  (N  .3 

).  2)  THEN 

WRITE  (*,*)  ' 

Enter  title  of  file  containing  parameters' 

READ  (M05) 

TTTL 

OPEN(l,FILE^l 

TTL) 

READ(1,105)  TITLE 

READ(1,110)  NPTS 

READ(1,110)  NRT 

READ(1,110)  Kd 

READ (1,110)  * 

[ 

READ (1,110)  DELTAY 

READ (1,110)  NSTRTPT 

READ (1,110)  N 

CAUS 

CLOSE(l) 

GO  TO  85 

ELSE 

GO  TO  16 

ENDIF 

WRITE  (*,*)  * 

• 

NUFILE=.TRUE. 

IF  (.NOT.  DSET)  NSTRTPT=1 

WRITE  (*,*)  'Enter  title  of  data  file  to  be  read' 

READ  (*,105)  TITLE 

OPEN  ( 1 ,  FILE=nTLE) 

READ(1,105)  HEADER 

READ  (1,110)  NPTS 

IF  (NPTS  .GT.  1024)  THEN 

WRITE  (*,*)  'Nianber  of  points  in  data  file  exceeds  the  dimension' 

WRITE  (*,*)  'of  the  array  used  in  the  program  to  store  the  file' 
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STOP 

ENDIF 

CLOSE(l) 

U  (DSET)  THEN 

IF  (NSTRTPT+(Rd+M-l)*DELTAY  .LE.  NPTS)  GO  TO  85 

ENDIF 

3  IF  (nufile)  THEN 

WRITE  (*,*)  'Enter  Rd,  >=  the  estimated  order  of  the  system  ' 

READ  (*,*)  Kd 

IF  (Rd  .GT.  69)  THEN 

WRITE  (*,*)  'Rd  must  be  less  than  70,  or  dimension  statements' 

WRITE  (*,*)  'in  this  program  must  changed  by  the  user' 

GO  TO  3 

ELSFJF  (Rd  .LT.  2)  THEN 

WRITE  (*,*)  'Rd  must  be  at  least  2' 

GO  TO  3 

ENDIF 
IF  (2*Rd  .GT.  NPTS)  THEN 

WRITE  (*,*)  'Rd  must  be  less  than  or  equal  to  \NPTS/2 
GOTO  3 

ELSEIF  (2*Rd  .EQ.  NPTS)  THEN 
WRITE  (*,*)  'Rd  equals' fRd 
WRITE  (*,*)  'M  must  be',Rd 
M=Rd 

WRITE  (*,*)  'since  there  are  a  total  of, NPTS 
WRITE  (*,*)  'points  in  ', TITLE 
GO  TO  45 
ENDIF 
GOTO  4 
ELSEIF  (DSET)  THEN 

N=M 
20     IF  (NSTRTPT+(N+M-1)*DELTAY  .LE.  NPTS)  THEN 

WRITE  (*,*)  'Given  the  other  parameters  chosen  thus  far,' 
25     WRITE  (*,*)  'Rd  may  range  from  \NRT 

WRITE  (*,*)  '  to',N 

WRITE  (*,*)  'Enter  Rd' 

READ  (*,*)  Rd 

IF  (Rd  .GE.  NRT  .AND.  Rd  .LE.  N)  GO  TO  85 

GOTO  25 

ELSE 

N=N-1 

GOTO  20 

ENDIF 
ENDIF 

4  IF  (NUFILE)  THEN 

WRITE  (*,*)  'Enter  M,  the  row  dimension  of  the  data  matrix' 
IF  (.NOT.  DSET  .AND.  LONG)  THEN 
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WRITE  (*,*)  '  ' 

WRITE  (*,*)  'Note:  Kd+M  points  in  ', title 
WRITE  (*,*)  '     will  be  processed  ' 
WRITE  (*,*)  '  * 
ENDIF 
30   WRITE  (*,*)  *M  may  range  from', Kd 
IF  (NPTS-Kd  .GT.  69)  THEN 
WRITE  (*,*)  '  to       69' 

ELSE 

WRITE  (*,*)  '  to* , NPTS-Kd 

ENDIF 

READ  (*,*)  M 
IF  (M  .GT.  69)  THEN 

WRITE  (*,*)  'M  must  also  be  less  than  70' 
GO  TO  30 

ELSEIF  (M  .LT.  Kd)  THEN 

WRITE  (*,*)  'M  must  be  greater  than  or  equal  to  Kd,  Kd=  \Kd 
GOTO  30 

ELSEIF  (Kd+tt  .GT.  NPTS)  THEN 

WRITE  (*,*)  'Kd4M  must  be  less  than  or  equal  to* ,NPTSf ' , ' 
WRITE  (*,*)  'the  number  of  data  points  in*, TITLE 
WRITE  (*,*)  '  ' 
GO  TO  30 
ENDIF 
ELSE 
N=Kd 
35    IF  (NSTRTPT+(Kd4N-l)*DELTAY  .LE.  NPTS)  THEN 

N=N+1 

GOTO  35 

ELSE 

N=N-1 

ENDIF 

IF  (N  .EQ.  Kd)  THEN 

WRITE  (*,*)  'M  must  equal', Kd 

M=Kd 

GOTO  85 

ENDIF 

IF  (N  .GT.  69)  N=69 
40   WRITE  (*,*)  'M  may  range  from', Kd 
WRITE  (*,*)  '  to\N 

WRITE  (*,*)  'Enter  M' 
READ  (*,*)  M 

IF  (M  .GE.  Kd  .AND.  M  .LE.  N)  GO  TO  85 
GO  TO  40 
ENDIF 

45    IF  (.NOT.  NUFILE)  GO  TO  85 

5    N=l 

50    IF  (NSTRTPT4N*(Kd-ftf-l)  .LE.  NPTS)  THEN 
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55 


N=N+1 

\ 

GOTO  50 

ELSE 

N=N-1 

ENDIF 

IF  (N  .EQ.  1)  THEN 

WRITE  (*,*) 

'Given  the  other  parameters  chosen  thus  far, ' 

WRITE  (*,*) 

'Spacing  can  only  be  1' 

DELTAY=1 

IF  (NUFILE) 

THEN 

GOTO  60 

ELSE 

GO  TO  85 

ENDIF 

ENDIF 

IF  (.NOT.  DSET  .AND.  LONG)  THEN 

WRITE  (*,*) 

'Enter  spacing  between  the  ' ,Kd+tf 

WRITE  (*,*) 

'data  points  of  ', TITLE 

WRITE  (*,*) 

'to  be  processed  ' 

WRITE  (*,*) 

■  i 

WRITE  (*,*) 

'If,  for  example,  one  is 

chosen,  then  \Kd4M 

WRITE  (*,*) 

'consecutive  points  in  ' 

, TITLE 

WRITE  (*,*) 

'will  be  processed  ' 

WRITE  (*,*) 

i  i 

ENDIF 

WRITE  (*,*) 

'Spacing  may  range  from 

1  ' 

WRITE  (*,*) 

to' 

,N 

READ  (*,*)  DELTAY 

IF  (DELTAY  . 

GE.  1  .AND.  DELTAY  .LE. 

N)  THEN 

IF  (NUFILE) 

THEN 

GOTO  60 

ELSE 

GOTO  85 

ENDIF 

ELSE 

GO  TO  55 

ENDIF 

60   WRITE  (*,*)  'Do  you  wish  to  adjust  eigenvalues?  (y/n) ' 

READ  (*,120)  YN 

IF  (YN  .EQ.  'N'  .OR.  YN  .EQ.  V)  THEN 

IF  (NUFILE)  GO  TO  6 

GOTO  85 

ENDIF 

IF  (YN  .NE. 
2    WRITE  (*,*) 

READ  (M20)  DC 

IF  (DC  .EQ.  'D'  .OR.  DC  .EQ.  'd')  GO  TO  65 

IF  (DC  .NE.  'C  .AND.  DC  .NE.  'c')  GO  TO  2 

WRITE  (*,*)  'Enter  estimate  of  the  actual  order  of  the  system' 


'Y'  .AND.  YN  .NE.  'y')  GO  TO  60 

'Discard  or  compensate  eigenvalues?  (d/c) ' 
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WRITE  (*,*)  '  ■ 
IF  (LONG)  THEN 

WRITE  (*,*)  'This  estimate  will  be  used  to  determine  the  ' 
WRITE  (*,*)  'number  of  eigenvalues  compensated  or  discarded  ' 
ENDIF 
65   WRITE  (*,*)  'the  estimate  may  range  from  2' 

WRITE  (*,*)  '  to',Kd-l 

READ  (*,*)  NRT 

IF  (NRT  .GT.  Kd  .OR.  NRT  .LT.  2)  THEN 
GOTO  65 

ELSEIF  (.NOT.  NUFILE)  THEN 
GOTO  85 
ENDIF 


.US.  NPTS)  THEN 


6    NSTRTPT=1 

70    IF  (NSTRTPT+(Kd-ftf-l)*DELTAY 

NSTRTPWISTRTPT+1 

GO  TO  70 

ELSE 

NSTRTPMISTRTPT-1 

ENDIF 

IF  (NSTRTPT  .EQ.  1)  THEN 

WRITE  (*,*)  'Given  the  other  parameters  chosen  thus  far,' 

WRITE  (*,*)  'the  starting  point  for  processing  the  data' 

WRITE  (*,*)  'must  be  the  first  point  in  the  data  file' 

GO  TO  85 

ENDIF 

WRITE  (*,*)  'Enter  desired  starting  point  in  data  file* 

IF  (.NOT.  DSET  .AND.  LONG)  THEN 

WRITE  (*,*)  '1  indicates  the  first  point  in  the  data  file 

ENDIF 

WRITE  (*,*)  '  ' 

WRITE  (*,*) 
75   WRITE  (*,*) 

WRITE  (*,*)  ' 

READ  (*,*)  N 

IF  (N  .GE.  1  .AND 

NSTRTPT=*J 

ELSE 

WRITE  (*,*)  'Enter  starting  point  again' 

WRITE  (*,*)  '  ' 

GOTO  75 

ENDIF 

IF  (.NOT.  NUFILE)  GO  TO  85 


'Given  the  other  parameters  chosen  thus  far, ' 
'the  starting  point  may  range  from       1' 

to', NSTRTPT 


N  .LE.  NSTRTPT)  THEN 


WRITE 

(* 

,*)  ' 

WRITE 

(* 

*)  ' 

WRITE 

(* 

*)  ' 

WRITE 

(*■ 

*)  ' 

WRITE 

(* 

*)  ' 

Do  you  want  the  data  matrix  arrangement  to  be' 

1.  Causal' 

2.  Non-causal' 
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80   WRITE  (*,*)  'Please  enter  1  or  2  ' 

READ  (*,*)  NCAUS 

IF  (NCAUS  .EQ.  1)  THEN 

CAUSAL=.TRUE. 

ELSEIF  (NCAUS  .EQ.  2)  THEN 

CAUSAL=.FALSE. 

ELSE 

GOTO  80 

ENDIF 

GOTO  85 
9    WRITE  (*,*)  'Enter  title  of  file  to  contain  parameters' 

READ  (*,105)  TETL 

OPEN(l,FILE==TnL) 

WRITE(1,105)  TITLE 

WRITE (1,110)  NPTS 

WRITE (1,110)  NRT 

WRITE(1,110)  Kd 

WRITE  (1,110)  M 

WRITE (1,110)  DELTAY 

WRITE  (1,110)  NSTRTPT 

WRITE  (1,110)  NCAUS 

CLOSE(l) 

IF  (DSET)  GO  TO  85 

12   IF  (DSET)  THEN 
CLOSE (2) 
CLOSEO) 

CALL  SUBPLT(NOVERLAY) 
ENDIF 


85   DSET=.TRUE. 

NUFILE=.  FALSE. 

WRITE(*,*)  '  ' 

WRITE(*,*)  '1.  Data  file  to  be  processed  ',T 

+ITLE 

WRTTE(*,*)  '  Number  of  data  points  in  data  file       \NPTS 

WRTTE(*,*)  '2.  Estimated  order  of  the  system  ' ,NRT 

WRITE(*,*)  '3.  Kd,  the  number  of  columns  in  the  data  matrix', Kd 

WRTTE(*,*)  '4.  M,  the  number  of  rows   in  the  data  matrix', M 

WRITE (*,*)  '5.  Spacing  between  data  points  being  processed  ', DELTA 

+Y 

WRITE(*,*)  '6.  First  point  in  the  data  file  to  be  processed', NSTRT 

+PT 

WRITE (*,*)  '   Last  point  in  the  data  file  to  be  processed', NSTRT 


IF 


+PT+Kd+M-1 

IF  (NCAUS  .BQ.  1)  THEN 

WRITE (*,*)  '7.  Data  matrix  arrangement  for  processing 
■KISAL   ' 

ELSE 


CA 
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WRITE(*,*)  '7.  Data  matrix  arrangement  for  processing     NON-CA 
-KJSAL   ' 

ENDIF 

WRITE(*,*)  "  ' 

WRITE(*,*)  '8.  Begin  processing  using  above  settings' 

WRTTE(*,*)  '9.  Store  parameters  1-7  in  a  file' 

WRITE (*,*)  '10.  Retrieve  parameters  1-7  from  a  previously  created 
+file' 

WRITE(*,*)  '11.  Reset  overlays' 

WRTTE(*,*)  '12.  Re-plot  overlays' 

WRTTE(*,*)  '13.  End  this  session  of  Kumaresan-Tufts  signal  process 
+ing' 

WRITE(*,*)  '  ' 

WRITE (*,*)  'Enter  an  integer  from  1  to  12  to  make  changes  as  often 
+  as  you  desire' 

90    READ  (*,*)  NMENU 

IF  (NMENU  .LT.  1  .OR.  NMENU  .GT.  13)  THEN 

WRITE(*,*)  'Enter  an  integer  from  1  to  13' 

GO  TO  90 

ENDIF 

GO  TO  (1,2,3,4,5,6,7,8,9,10,11,12,13), NMENU 

8    OPEN(l,FILE=nni:) 

READ (1,105)  HEADER 

READ (1,110)  NPTS 

READ (1,115)  XQ 

READ (1,115)  XQ 

DO  95  1=1, NPTS 

READ(1,115)  D(I) 
95    CONTINUE 

CLOSE(l) 

KdPLT=Kd 

WRITE (*,*)  'Enter  title  of  file  to  contain  real  part  of  poles' 

READ (*,  105)  TnUR 

OPEN  (2 , f ile=TTTLER) 

WRITE(*,*) 'Enter  title  of  file  to  contain  imaginary  part  of  poles' 
READ(*,105)  Trnn 
0PEN(3,file=TnTJEI) 
WRITE (10, 100)    (KdPLT) 
WRTTE(10,105)  TTTLER 
WRITE(10,105)  TTTLEI 
100      FORMAT (12) 

MN=MAX(M,Kd) 

105      FORMAT  (A) 
110      FORMAT (15) 
115      FORMAT  (FJ.2. 6) 
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120   FORMAT  (Al) 
C    Form  data  matrix 
DO  125  I=l,Kd+M 

dy(i)=d(  (1-1)  *deltay*nstrtpt) 
125    continue 

130   DO  140  1=1, M 

DO  135  J=l,Kd 

A(I,J)=Dy(I+J) 
135   CONTINUE 
140   CONTINUE 

B(l)=Dy(l) 
DO  145  1=2, M 
B(I)=A(I-1,1) 
145   CONTINUE 

C    Begin  singular  value  decomposition 

CALL  SVD(MACHEP,M,Kd,MN,A,W,MATU,U,MATV,V,IERR,RVl) 

C    Errors  in  SVD? 

IF  (IERR  .GT.  0.0)  THEN 

WRITE  (*,*)  'Error  in  singular  value  number  ', IERR, STOP 

ENDIF 

IF  (YN  .EQ.  'N')  GO  TO  190 

DO  150  1=1, Kd 
XP(I)=0.0 
150   CONTINUE 

C    Discard  or  compensate  eigenvalues 
C    Order  singular  values 

XP(1)=W(1) 

DO  165  1=2, Kd 

DO  160  J=1,I 

IF  (W(I)  .GT.  XP(J))  THEN 

DO  155  K=I+1,J,-1 
155   XP(K)=XP(K-1) 

XP(J)=W(I) 

GO  TO  165 

ENDIF 
160   CONTINUE 

XP(I+1)=W(I) 
165   CONTINUE 

C    XP{  )  now  contains  ordered  singular  values-XP(l)  is  the  largest 
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c 

Discard  eigenvalues 

IF  (DC  .ECj.  'D')  THEN 

DO  170  J=NRT+l,Kd 

170 

W(J)=(0.0) 

ELSE 

C 

Compensate  eigenvalues 

AVG=0.0 

DO  175  J=NRT+l,Kd 

AVG=AVG+XP(J)**2 

175 

CONTINUE 

IF  (Rd  .GT.  NRT)  AVG=AVG/DBLE (FLOAT (Kd-NRT) 

DO  185  J=l,Kd 

DO  180  K=l,Kd 

IF  (  W(J)  .ECj.  XP(K)  )  THEN 

IF  (  K  .GT.  NRT  )  THEN 

W(J)=0.0 

ELSE 

V(J)=DSOj*T(DABS(  W(J)*W(J)-AVG)) 

ENDIF 

GOTO  185 

ENDIF 

180 

CONTINUE 

185 

CONTINUE 

ENDIF 

190 

DO  200  1=1, M 

DO  195  J=1,M 

UT(I,J)=(U(J,I)) 

195 

CONTINUE 

200 

CONTINUE 

c 

Form  SIGMA+  (KflxM) 

DO  210  1=1 ,Kd 

DO  205  J=1,M 

SIGMA (I, J) =0.0 

IF  (I  .EQ.  J  .AND.  W(J)  .NE.  0.0)  THEN 

SIGMA(I,J)=1.0D0/V(J) 

ELSE 

SIGMA(I,J)=O.0d0 

ENDIF 

205 

CONTINUE 

210 

CONTINUE 

Form  SIGMA  (MxKd) 

DO  220  I=1,M 

DO  215  J=l,Kd 

SIG(I,J)=0.0 

IF  (I  .EQ.  J)  SIG(I,J)=W(J) 
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215      camMJE 

220        CONTINUE 

C         V=KdxKd,SI<m+=KdxM,VS=KdxM 

CAIL  MXMUL(V,SIGMA,Kd,Kd,M,VS) 

C         TO=KdxM,trT=MxM,AINV=KdxM 

CALL  MJMJL(VS,UT,Kd,M,M,AINV) 

C    Calculate  matrix  multiplication  of  AINV  x  B,  where 
C    AINV=KdxM,B=Mxl,XP=Kdxl 

CALL  MXMUL(AINV,B,Kd,M,L,XP) 

C    Calculate  autoregressive  coefficients  from  prediction  coefficients 

IF  (XP(Kd)  .ECj.  0.0)  THEN 

WRITE  (*,*)  'ERROR,  avoiding  division  by  zero' 

STOP 

ELSE 

B(Kd)=1.0d0/XP(Kd) 

ENDIF 

DO  225  1=2, Kd 

B (1-1 ) =-B (Kd) *XP (Kd-I+1 ) 
225   CONTINUE 

DO  230  1=1,  Kd 
X (I) =-B (Kd-I+1) 

IF  (NCAUS  .BQ.  1)  X (I) =-XP (Kd-I+1) 
230   CONTINUE 

X(Kd+l)=1.0 

C    Compute  the  roots  of  the  polynomial  in  z 

CALL  P0LRT(X/C0F,KD,RCOTR,RCOTI,IER) 

IF  (IER  .NE.  0)  WRITE  (*,*)  'ERROR  with  POLRT,  IER=' ,IER,STOP 

DO  235  1=1, Kd 
WRTTE(2,115)  ROOTR(I) 
WRTTE(3,115)  RCOn(I) 
S (I)=DCMPLX(RCCTR(I) ,RCOTI (I) ) 
235   CONTINUE 

MAGPOL=0 
DO  240  1=1, Kd 

IF  (CDABS(S(D)  .GE.  l.OdO)  MAGP0L=MAGP0L+1 
240   CONTINUE 

WRTTE(*,*)  '#  of  poles  with  magnitude  <=  l',Kd-MAGPOL 
WRITE  (*,*)  'HIT  ANY  KEY  TO  CONTINUE' 
READ  (*,105)  HEADER 
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C    Plot  poles 

N0VERIAY=NOVERLAY+1 

CLOSE  (2) 

CLOSED) 

CALL  SUBPLT(NOVERLAY) 

J=0 
K=0 

DO  245  1=1, Kd 

IF  (CDABS(S(I))  .LT.  1.0)  THEN 

J=J+1 

K=K+1 

WRITE  (*,*)  S(I),CDABS(S(I)) 

ENDIF 

IF  (J  .BQ.  20)  THEN 

WRITE  (*,*)  'Enter  any  key  to  continue' 

READ  (*,105)  HEADER 

J=0 

ENDIF 

245  continue 

WRTTE(*f*)  'Poles  with  magnitude  less  than  one:  \K 

GO  TO  85 
13   STOP 
END 
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APPENDIX  B:   THE  CADZOW-SOLOMON  POLE  EXTRACTION  ALGORITHM 

The  following  program  implements  the  Cadzow-Solomon 
algorithm  as  described  in  Chapter  2  of  this  thesis.  The 
program  is  written  in  Fortran  77.  The  SVD  and  root-finding 
subroutines  called  by  this  program  are  found  in  the  EISPACK 
library  [18] .  The  SVD  subroutine  is  a  translation  from  ALGOL 
as  given  in  [19]  .  The  matrix  multiplication  and  graphics 
subroutines,  also  called  by  this  program,  are  found  in 
Appendix  C  and  D  respectively. 


SLARGE 

INTEGER  IERR,Kd,Kn,M,m,MAGPOL,NSTRTPT,DELTAY 

INFECT*  IER,NCAUS,NMBNU,INSTRTPT 

HfTEGER*2  KdPLT 

REAL*8  A(70,70) ,W(70) ,U(70,70) ,V(70,70) ,RV1(70) 

REAL*8  VS (70,70) ,UT(70,70) ,AINV(70,70) ,X(70) 

REAL*8  XP(70) ,B(70) ,SIGMA(70,70) ,SIG(70,70) 

REAL*8  OCF(70),ROOTR(70),ROOTI(70) 

REAL   MAG 

REAL*8  D(1024) ,AVG,MACHEP/1.0E-16/,Dy(14O) ,Dx(1024) 

00MPLEX*16  S(70) 

LOGICAL  MATU/.TOUE7,MATV/.TRUE./,CAUSAL/.TOUE./, LONG/. TRUE./ 

LOGICAL  DSET/. FALSE. /,NUFILE/. TRUE./ 

CHARACTER  TITIJE*16,HEADER*64,YNn,DC*l,TITLERn6,TITLEI*16 

CHARACTER  TITL*16,TrrLD*16 

C         Enter  parameters  for  processing 

14        IF  (DSET)  CLOSE  (10) 
NOVERLAY=0 
OPEN(10,FILE=,PLOT') 
IF  (DSET)  GO  TO  215 

WRITE  (*,*)   'Welcome  to  signal  processing  using  the' 
WRITE  (*,*)    'Cadzow-Solomon  method' 
WRITE  (*,*)    '    ' 
WRITE  (*,*)   'Do  you  want  ' 
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WRITE  (*,*) 
WRITE  (*,*) 
WRITE  (*,*) 
WRITE  (*,*) 
25  WRITE  (*,*) 
READ  (*,*)  N 
IF  (N  .EQ.  1)  THEN 
LCNG=.TRUE. 

ELSEIF  (N  .BQ.  2)  THEN 
LONG=.  FALSE. 
ELSE 
GOTO  25 
ENDIF 


1.  The  long  version  for  beginners' 

2.  The  short  version  for  pros' 

Please  enter  1  or  2  ' 


35 


13 


WRITE  (*,*) 
+r  processing 
WRITE  (*,*) 
(* 


{* 


WRITE 
WRITE 
WRITE  (* 
WRITE  (* 
WRITE  (* 


.*) 
.*) 
*) 
,*) 
*) 


Session  will  begin  with  entry  of  parameters  needed  fo 


Do  you  want  to  enter  parameters  from' 

1.  The  keyboard' 

2.  A  previously  created  file  of  parameters' 
i 

Please  enter  1  or  2  ' 


WRITE  (*,*) 

READ  (*,*)  N 

IF  (N  .EQ.  1)  THEN 

GO  TO  8 

ELSEIF  (N  .EQ.  2)  THEN 

WRITE  (*,*)  'Enter  title  of  file  containing  parameters' 

READ  (MOO)  TTTL 

OPEN(l,FILE=nTL) 

READ(1,100)  TITLE 

READ (1,110)  NPTS 

READ(1,110)  NRT 

READ (1,110)  Kd 

READ (1,110)  M 

READ (1,110)  DELTAY 

READ (1,110)  NSTRTPT 

READ (1,110)  NCAUS 

READ(1,100)  TITLD 

READ(1,110)  NDPTS 

READ (1,110)  Kn 

READ (1,110)  INSTRTPT 

CLOSE(l) 

GO  TO  215 

ELSE 

GO  TO  35 

ENDIF 

WRITE  (*,*)  '  ' 


8    WRITE  (*,*)  'Enter  title  of  file  containing  excitation  waveform' 
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READ  (\100)  TTTLD 

OPEN  (8 ,  FILE=nTLiD) 

READ(8,100)  HEADER 

READ(8,110)  N 

IF  (N  .GT.  1024)  THEN 

WRITE  (*,*)  'Number  of  points  in  data  file  exceeds  the  dimension' 

WRITE  (*,*)  "of  the  array  used  in  the  program  to  store  the  file' 

STOP 

ENDIF 

CLOSE(8) 

IF  ((N  .GE.  NDPTS)  .AND.  DSET)  THEN 

NDPTS=N 

GOTO  215 

ENDIF 

NDPTS=N 

9  WRITE  (*,*)  'Enter  estimated  order  of  waveform' 
IF  (DSET)  THEN 

MAXIMUM=NDPTS-M 

IF  (MAXIMUM  .GT.  M-Kd-1)  MAXIMUM=M-Kd-1 
IF  (MAXIMUM  .GT.  NDPTS-INSTRTPT-Kn-M+1)  THEN 
MAXIMUhH©PTS-INSTRTPT-Kn-M+l 
ENDIF 
ELSE 

MAXIMUM=66 
ENDIF 

IF  (MAXIMUM  .EQ.  1)  THEN 

WRITE  (*,*)  'The  estimated  order  of  the  waveform  can  only  be  1' 
IF  (DSET)  GO  TO  215 
GO  TO  10 
ELSE 

IF  (DSET)  THEN 

WRITE  (*,*)  'Given  the  other  parameters  chosen  thus  far,' 
ENDIF 
45   WRITE  (*,*)  'the  order  may  range  from  1' 

WRITE  (*,*)  '  to', MAXIMUM 

READ  (*,*)  Kn 

IF  (Kn  .GE.  1  .AND.  Kn  .IE.   MAXIMUM)  THEN 
IF  (DSET)  GO  TO  215 
GO  TO  10 
ENDIF 

WRITE  (*,*)  'Enter  estimated  order  again* 
WRITE  (*,*)  '  ' 
GO  TO  45 
ENDIF 
IF  (DSET)  GO  TO  215 

10  INSTRTPT=1 

55        IF  (mSTRTPT+Kn+M-l  .GT.  NDPTS)  THEN 
INSTRTFI-INSTRTPT-1 
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ELSE 

INSTRTPT=INSTRTPT+1 

GO  TO  55 

ENDIF 

MSTRT=INSTRTPT 

IF  (INSTRTPT  .EQ.  1)  THEN 
WRITE  (*,*)  'The  first  point  can  only  be  1' 
GOTO  215 
ELSE 

WRITE  (*,*)  'Enter  first  point  in  waveform  file  to  be  processed' 
65   WRITE  (*,*)  'Given  the  other  parameters  chosen  thus  far,' 
WRITE  (*,*)  'the  starting  point  may  range  from       1' 
WRITE  (*,*)  '  to\MSTRT 

READ  (*,*)  INSTRTPT 

IF  (INSTRTPT  .GE.  1  .AND.  INSTRTPT  .LE.  MSTRT)  THEN 
IF  (DSET)  GO  TO  215 
GO  TO  1 
ENDIF 

WRITE  (*,*)  'Enter  starting  point  again' 
WRITE  (*,*)  '  ' 
GOTO  65 
ENDIF 
IF  (DSET)  GO  TO  215 

1    IF  (.NOT.  DSET)  NUFILE=.TRUE. 
IF  (.NOT.  DSET)  NSTRTPT=1 

WRITE  (*,*)  'Enter  title  of  data  file  to  be  read' 
READ  (MOO)  TITLE 
OPEN (12 , FTLE=TrTLE) 
READ (12, 100)  HEADER 
READ (12, 110)  NPTS 
IF  (NPTS  .GT.  1024)  THEN 

WRITE  (*,*)  'Number  of  points  in  data  file  exceeds  the  dimension' 
WRITE  (*,*)  'of  the  array  used  in  the  program  to  store  the  file' 
STOP 
ENDIF 
CL0SE(12) 

IF  (NUFILE)  THEN 

GOTO  3 

ELSEU  (NSTRTPT+(Kd4M-l)*DELTAY  .LE.  NPTS)  THEN 

GO  TO  215 

ELSE 

GO  TO  6 

ENDIF 


IF  (NUFILE)  THEN 
MAXMJM=69-Kn-1 
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IF  (MAXIMUM  .GT.  NPTS-69)  MAXIMUM=NPTS-69 

MD*=2 

IF  (MIN  .EQ.  MAXIMUM)  THEN 

Kd=MIN 

WRITE  (*,*)  'Given  the  other  parameters  chosen  thus  far,' 

WRITE  (*,*)  'Kd  must  be  \MIN 

GOTO  4 

ENDIF 

WRITE  (*,*)  'Enter  Kd,  >=  the  estimated  order  of  the  system  ' 

WRITE  (*,*)  'Given  the  other  parameters  chosen  thus  far,' 
75   WRITE  (*,*)  'Kd  may  range  franMDN 

WRITE  (*,*)  '  to', MAXIMUM 

READ  (*,*)  Kd 

IF  (Kd  .GE.  MIN  .AND.  Kd  .LE.  MAXIMUM)  GO  TO  4 
GO  TO  75 

ELSFJF  (DSET)  THEN 
MAXIMUM=M-Kn-1 

IF  (MAXIMUM  .GT.  NPTS-M)  MAXIMUM=NPTS-M 
MIN=2 
N=MAXIMUM 
85    IF  (NSTRTPT+(N+M-1)*DELTAY  .LE.  NPTS)  THEN 
MAXIMUMS 

IF  (MIN  .EQ.  MAXIMUM)  THEN 
Kd=MIN 
GOTO  215 

ELSFJF  (MAXIMUM  .LT.  MIN)  THEN 
DELTAY=1 

IF  (1+(24M-1)*DELTAY  .LE.  NPTS)  THEN 
Kd=2 

GO  TO  135 
ENDIF 

WRITE  (*,*)  'Error.  Kd  must  be  less  than  2' 
Kd=2 

GO  TO  215 
ENDIF 

WRITE  (*,*)  'Given  the  other  parameters  chosen  thus  far,' 
95     WRITE  (*,*)  'Kd  may  range  from  \MIN 

WRITE  (*,*)  '  to', MAXIMUM 

WRITE  (*,*)  'Enter  Kd' 

READ  (*,*)  Kd 

IF  (Kd  .GE.  MIN  .AND.  Kd  .LE.  MAXIMUM)  GO  TO  215 

GO  TO  95 
ELSE 

N=N-1 

GOTO  85 
ENDIF 
ENDIF 
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C    Determine  M 

4    IF  (NUFILE)  THEN 

WRITE  (*,*)  'Enter  M,  the  row  dimension  of  the  data  matrix1 
IF  (.NOT.  DSET  .AND.  LONG)  THEN 
WRITE  (*,*)  '  ' 

WRITE  (*,*)  'Note:  Kd+M  points  in  ', title 
WRITE  (*,*)  '     will  be  processed  ' 
WRITE  (*,*)  '  ' 
ENDIF 
105   WRITE  (*,*)  'M  may  range  fran',Kd 
IF  (NPTS-Kd  .GT.  69)  THEN 
WRITE  (*,*)  '  to       69' 

ELSE 

WRITE  (*,*)  '  to', NPTS-Kd 

ENDU 

READ  (*,*)  M 
IF  (M  .(7T.  69)  THEN 

WRITE  (*,*)  'M  must  also  be  less  than  70' 
GOTO  105 

ELSEIF  (M  .LT.  Kd)  THEN 

WRITE  (*,*)   'M  must  be  greater  than  or  equal  to  Kd,  Kd=  ',Kd 
GOTO  105 

ELSEIF  (Kd+M  .GT.  NPTS)  THEN 

WRITE  (*,*)   'Kd-tfl  must  be  less  than  or  equal  to' ,NPTS, ' , ' 
WRITE  (*,*)   'the  number  of  data  points  in',TTTLE 
WRITE  (*,*)    '   ' 
GOTO  105 
ENDIF 
C         Begin  part  for  data  already  set 
ELSE 
N=Kd 
115  IF  (NSTRTPT+(Kd-fN-l)*DELTAY  .LE.  NPTS)  THEN 

N=«+l 

GO  TO  115 

ELSE 

N=N-1 

ENDIF 

IF  (N  .EQ.  Kd)  THEN 

WRITE  (*,*)    'M  must  equal', Kd 

M=Kd 

GOTO  215 

ENDIF 

MAXDUHt 

IF  (MAXIMUM  .GT.  69)  MAXHHJM=69 

IF  (Kd+Kn+1  .EQ.  MAXIMUM)  THEN 

M=Kd4Kn+l 

GO  TO  215 

ELSEIF  (Kd+Kn+1  .GT.  MAXIMUM)  THEN 

WRITE  (*,*)    'Kd  must  be  reduced' 


130 


00  TO  3 

ELSE 

MDJ=Kd+Kn+l 

ENDIF 
IF  (MIN  .LT.  Kn+Kd+1)  MIN=Kn+Kd+l 
125   WRITE  (*,*)  'M  may  range  from', MIN 

WRITE  (*,*)  '  to', MAXIMUM 

WRITE  (*,*)  'Enter  M* 

READ  (*,*)  M 

IF  (M  .GE.  MIN  .AND.  M  .LE.  MAXIMUM)  GO  TO  215 

GO  TO  125 

ENDIF 


c    Determine  DELTAY 

135   IF  (.NOT.  NUFTLE)  GO  TO  215 

5    N=l 

145   IF  (NSTRTPT+NMftMKL)  .LE.  NPTS)  THEN 

N=N+1 

GO  TO  145 

ELSE 

N=fl-1 

ENDIF 

IF  (N  .EQ.  1)  THEN 

WRITE  (*,*)  'Given  the  other  parameters  chosen  thus  far, 

WRITE  (*,*)  'Spacing  can  only  be  1' 

DELTAY=1 
IF  (NUFILE)  THEN 
GO  TO  165 
ELSE 

GO  TO  215 
ENDIF 

ENDIF 

IF  (.NOT 

WRITE  (* 

WRITE  (* 

WRITE  (* 

WRITE  (* 

WRITE  (* 

WRITE  (* 

WRITE  (* 

WRITE 

ELSE 

WRITE 

WRITE  (* 

ENDIF 
155  WRITE  (*,*)  'Spacing  may  range  from 

WRITE  (*,*)  '  to',N 

READ  (*,*)  DELTAY 

IF  (DELTAY  .GE.  1  .AND.  DELTAY  .LE.  N)  THEN 
IF  (NUFILE)  THEN 


[* 


0 


DSET  .AND.  LONG)  THEN 


0 


Enter  spacing  between  the  ',Kd+M 

data  points  of  ' , TITLE 

to  be  processed  ' 
i 

If,  for  example,  one  is  chosen,  then  ',Kd+M 
consecutive  points  in  ' , TITLE 
will  be  processed  ' 


Enter  spacing 
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GOTO  165 
ELSE 

GOTO  215 
ENDIF 

T1.SE 

GOTO  155 
ENDIF 

165  WRITE  (*,*)  'Do  you  wish  to  adjust  eigenvalues?  (y/n) ' 

READ  (M50)  YN 

IF  (YN  .EQ.  'N'  .OR.  YN  .EQ.  'n')  THEN 

IF  (NUFILE)  GO  TO  6 

GO  TO  215 

ENDIF 

IF  (YN  .NE.  'Y'  .AND.  YN  .NE.  'y')  GO  TO  165 
2    WRITE  (*,*)  'Discard  or  compensate  eigenvalues?  (d/c) ' 

READ  (\150)  DC 

IF  (DC  .EQ.  'D'  .OR.  DC  .EQ.  'd')  THEN 

NRT=Kd 

GOTO  175 

ENDIF 

IF  (DC  .NE.  'C  .AND.  DC  .NE.  'c')  GO  TO  2 

WRITE  (*,*)  'Enter  estimate  of  the  actual  order  of  the  system' 

WRITE  (*,*)  '  ' 

IF  (LONG)  THEN 

WRITE  (*,*)  'This  estimate  will  be  used  to  determine  the  ' 

WRITE  (*,*)  'number  of  eigenvalues  compensated  or  discarded  ' 

ENDIF 
175   WRITE  (*,*)  'the  estimate  may  range  from  2' 

WRITE  (*,*)  '  to\Kd+Kn+l 

READ  (*,*)  NRT 

IF  (NRT  .GT.  Kd+Kn+1  .OR.  NRT  .LT.  2)  THEN 

GO  TO  175 

ELSEEF  (.NOT.  NUFILE)  THEN 

GO  TO  215 

ENDIF 

6    NSTRTPT=1 

185   IF  (NSTRTPT+(KcWI-l)*DELTAY  .LE.  NPTS)  THEN 

NSTRTPT=fKTRTPT+l 

GO  TO  185 

ELSE 

NSTRTPWJSTRTPT-1 

ENDIF 

IF  (NSTRTPT  .EQ.  1)  THEN 

WRITE  (*,*)  'Given  the  other  parameters  chosen  thus  far,' 

WRITE  (*,*)  'the  starting  point  for  processing  the  data' 

WRITE  (*,*)  'must  be  the  first  point  in  the  data  file* 

GO  TO  215 

ENDIF 
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WRITE  (*,*)  'Enter  desired  starting  point  in  data  file' 
IF  (.NOT.  DSET  .AND.  LONG)  THEN 

WRITE  (*,*)  '1  indicates  the  first  point  in  the  data  file 
ENDIF 

WRITE  (*,*)  '  ' 

WRITE  (*,*)  'Given  the  other  parameters  chosen  thus  far,' 
195  WRITE  (*,*)  'the  starting  point  may  range  from       1' 
WRITE  (*,*)  '  to',NSTRTPT 

READ  (*,*)  N 

IF  (N  .GE.  1  .AND.  N  .LE.  NSTRTPT)  THEN 
NSTRTPT=N 
ELSE 

WRITE  (*,*)  'Enter  starting  point  again' 
WRITE  (*,*)  '  ' 
GOTO  195 
ENDIF 
IF  (.NOT.  NUFILE)  GO  TO  215 


205 


IF  (DSET)  THEN 

IF  (NCAUS  .EQ.  1)  THEN 

NCAUS=2 

GOTO  215 

ELSE 

NCAUS=1 

GOTO  215 

ENDIF 

ENDIF 

WRITE  (* 

WRITE  (* 

WRITE  (* 

WRITE  (* 

WRITE  (* 

WRITE  (* 

READ 


Do  you  want  the  data  matrix  arrangement  to  be' 

1.  Causal' 

2.  Non-causal' 


Please  enter  1  or  2 
*,*)  NCAUS 


IF  (NCAUS  .EQ. 
CAUSAL=.TRUE. 
ELSETF  (NCAUS 
CAUSALf.  FALSE. 
ELSE 

GOTO  205 
ENDIF 
GOTO  215 


1)  THEN 


.EQ.  2)  THEN 


12   WRITE  (*,*)  'Enter  title  of  file  to  contain  parameters' 
READ  (MOO)  TTTL 
OPEN  ( 1 ,  FTLE==TTrL) 
WRITE(1,100)  TITLE 
WRTrE(l,110)  NPTS 
WRITE (1,110)  NRT 
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WRITE  (1,110) 

write(i,iio) 
VRrre(i,iio) 

WRITE  (1,110) 
WRITE (1,110) 

WRrre(i,ioo) 
write  (i,no) 

WRITE  (1,110) 
WRITE  (1,110) 
CLOSE(l) 
IF  (DSET)  GO 


Kd 

M 

DELTAY 

NSTRTPT 

NCAUS 

TTTLD 

NDPTS 

Kn 

BISTRTPT 

TO  215 


15   IF  (DSET)  THEN 
CL0SE(2) 
CL0SE(3) 

CALL  SUBPLT(NOVERLAY) 
ENDIF 


215   DSET=.TRUE. 

NUFILE=. FALSE. 
WRTTE(*,*)  '  ' 
WRITE(*,*) 

+ITLE 
WRITER,*) 
WRITE(*,*) 
WRTIE(*,*) 
WRTTE(*,*)  '4 
WRITE**,*)  '5 

+Y 
WRTIE(\*) 

+PT 
WRTTE(*,*) 

+PT+Kd+M-1 
IF  (NCAUS 
WRITER,*) 

■KJSAL   ' 
ELSE 
WRTTE(*,*)  '7 

-KJSAL   ' 
ENDIF 
WRITE(*,*) 


'1.  Data  file  to  be  processed 


',T 


'2. 
•3. 


•6. 


.EQ. 
'7. 


WRTTE(*,*) 
+ITLD 
WRTTE(*,*) 
WRITE(*,*) 

WRnE(*,*) 
WRTTE(*,*) 
+TPT 


Number  of  data  points  in  data  file       ' ,NPTS 
Estimated  order  of  the  system  ' ,NRT 

Rd,  the  number  of  columns  in  the  data  matrix', Kd 
M,  the  number  of  rows   in  the  data  matrix' ,M 
Spacing  between  data  points  being  processed  ' , DELTA 

First  point  in  the  data  file  to  be  processed', NSTRT 

Last  point  in  the  data  file  to  be  processed', NSTRT 


1)  THEN 

Data  matrix  arrangement  for  processing 


Data  matrix  arrangement  for  processing 


CA 


NON-CA 


8.  File  containing  excitation  waveform 

Number  of  data  points  in  above  file 

9.  Estimated  order  of  the  waveform 

10.  First  point  in  the  file  to  be   ' 
input  into  the  data  matrix 


',T 

', NDPTS 
',Kn 


MNSTR 
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WRTTE(*,*)  '  ' 

WRITE(*,*)  '11.  Begin  processing  using  above  settings' 

WRTTE(*,*)  '12.  Store  parameters  1-10  in  a  file' 

WRITE (*,*)  '13.  Retrieve  parameters  1-10  from  a  previously  created 
+  file' 

WRITE(*,*)  '14.  Reset  overlays' 

WRITE (*,*)  '15.  Re-plot  overlays' 

WRITE (*,*)  '16.  End  this  session  of  Cadzow-Solomon  signal  processi 
+ng' 

WRTTE(*,*)  '  ' 

WRITE (*,*)  'Enter  an  integer  from  1  to  16  to  make  changes  as  often 
+  as  you  desire' 
225   READ  (*,*)  NMENU 

IF  (NMENU  .LT.  1  .OR.  NMENU  .GT.  16)  THEN 

WRITE (*,*)    'Enter  an  integer  from  1  to  16' 

GO  TO  225 

ENDIF 

GO  TO  (1,2, 3, 4, 5, 6, 7, 8, 9,10, 11, 12, 13, 14, 15, 16), NMENU 

11        0PEN(12,FILE=TrTI£) 

READ(12,100)  HEADER 

READ(12,110)  NPTS 

READ (12, 120)  XQ 

READ (12, 120)  XQ 

DO  235  1=1, NPTS 

READ(12,120)  D(I) 
235      CONTINUE 

CljOSE(12) 

0PEN(8,FILE=TTTLD) 
READ (8, 100)  HEADER 
READ(8,110)  NDPTS 
READ (8, 120)  XQ 
READ (8, 120)  XQ 
DO  245  1=1, NDPTS 
READ(8,120)  Dx(I) 
245   CONTINUE 
CL0SE(8) 

KdPLT=Kd 

WRITE (*,*)  'enter  title  of  file  to  contain  real  part  of  poles' 

READ(MOO)  TITLER 

OPEN(2,FTLE=nTLER) 


WRTTE(*,*) 'enter  title  of  file  to  contain  imaginary  part  of  poles' 

READ(MOO)  TITLEI 

0PEN(3,nLE=TTTLEI) 
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WRITE (10, 130)    (KdPLT) 
WRITE (10, 100)  TTTLER 

write  do,  ioo)  tttlei 

130   FORMAT  (12) 

MN=MAX(M,Kd+Kn+l) 

100  FORMAT  (A) 

110  FORMAT (15) 

120  FORMAT (E12. 6) 

150  FORMAT  (A) 


DO  255  I=l,Kd4M 
Dy (I) =D ( (1-1 ) *DELTAY4flSTRTPT) 
255   CONTINUE 

265   DO  285  1=1, M 

DO  275  J=l,Kd+Kn+l 

A(I,J)=Dy(I+J) 

IF  (J  .GE.  Kd+1)  A(I,J)=Dx(I-KJ+INSTRTPT-2-Kd) 


275 

CONTINUE 

285 

CONTINUE 

B(l)=Dy(l) 

DO  295  I=2,M 

B(I)=A(I-1,1) 

295 

OCNTINUE 

N=Kd+Kn+l 

C    Begin  singular  value  decomposition 

CALL  SVD(MACHEP,M,N,MN,A,W,MATU,U,MATV,V,IERR,RV1) 

C    Errors  in  SVD? 

IF  (IERR  .GT.  0.0)  THEN 

WRITE  (*,*)  'Error  in  singular  value  number  ', IERR, STOP 

ENDIF 

IF  (YN  .EQ.  'N')  GOTO  385 

DO  305  I=l,Kd+Kn+l 
XP(I)=0.0 
305   CONTINUE 

C    Discard  or  compensate  eigenvalues 
c    Order  singular  values 

XP(1)=W(1) 
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DO  335  I=2,Kd+Kn+l 

DO  325  J=1,I 

if  (W(I)  .GT.  XP(J))  THEN 

DO  315  K=I+1,J,-1 
315   XP(K)=XP(K-1) 

XP(j)=W(i) 

GOTO  335 

ENDIF 
325   CONTINUE 

XP(I+1)=V(I) 
335   CONTINUE 

C    XP(  )  now  contains  ordered  singular  values:  XP(1)  is  the  largest 

C    Discard  eigenvalues 

IF  (DC  .EC;.  'D')  THEN 

DO  345  J=NRT+1,  Kd+Kn+1 
345   W(J)=(0.0) 

ELSE 
C    Compensate  eigenvalues 

AVG=0.0 

DO  355  J=NRT+1, Kd+Kn+1 

AVG=AVG+XP(J)**2 
355   CONTINUE 

IF  (Kd+Kn+1  .GT.  NRT)  AVG=AVG/DBLE (FLOAT (Kd+Kn+1-NRT) ) 

DO  375  J=l, Kd+Kn+1 
DO  365  K=l, Kd+Kn+1 
IF  (  W(J)  .EQ.  XP(K)  )  THEN 

IF  (  K  .GT.  NRT  )  THEN 

V(J)=0.0 

ELSE 

V(J)=DSQRT(DABS(  V(J)*W(J)-AVG)) 

ENDIF 

GOTO  375 

ENDIF 
365   CONTINUE 
375   CONTINUE 
ENDIF 


385   DO  405  1=1, M 

DO  395  J=1,M 

UT(I,J)=(U(J,I)) 
395   CONTINUE 
405   CONTINUE 

C    Form  SIGMA+  (Kd+Kn+1  x  M) 
DO  425  1=1, Kd+Kn+1 
DO  415  J=1,M 
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SIGMA(I,J)=0.0 

IF  (I  .EQ.  J  .AND.  W(J)    .NE.  0.0)  THEN 

SIGMA(I,J)=1.0d0/W(J) 

ELSE 

SIGMA(I,J)=O.0D0 

ENDIF 
415        CONTINUE 
425        CONTINUE 

C         Form  SIGMA  (M  x  Kd-tfCn+1) 

DO  445  1=1, M 

DO  435  J=l,Kd+Kh+l 

SIG(I,J)=0.0 

IF  (I  .EQ.  J)  SIG(I,J)=V(J) 
435   CONTINUE 
445   CONTINUE 

C  V=Kd4Kn+lxKd+KiHl,SIGto^ 

CALL  MXMUL(V,SIGMA,K(Mn+l,Kd-H(n+l,M,VS) 

C    VS=Kd+Kn+lxM,UT=MxM,AINVN(d-H(n+lxM 
CALL  MXMUL(VS,UT,Kd4Kn+l,M,M,AINV) 

C    Calculate  matrix  multiplication  of  AINV  x  B,  where 
C    AINV=Kd+KrH-lxM,B=Mxl,XP=Kd+Kn+lxl 
CALL  MXMUL(AINV,B,Kd-H<n+l,M,L,XP) 

C    Compute  autoregressive  coefficients  from  prediction  coefficients 

IF  (XP(Kd)  .EQ.  0.0)  THEN 
WRITE  (*,*)  'ERROR,  avoiding  division  by  zero' 
STOP 
ELSE 

B(Kd)=1.0d0/XP(Kd) 
ENDIF 

DO  455  1=2, Kd 
B (1-1) =-B (Kd) *XP (Kd-i+1) 
455   CONTINUE 

DO  465  i=l,Kd 
X(I)=-B(Kd-I+l) 

IF  (NCAUS  .EQ.  1)  X(I)=-XP(Kd-I+l) 
465   CONTINUE 

X(Kd+l)=1.0 

C    Compute  the  roots  of  the  polynomial  in  z 

CALL  POLRT(X,COF,KD,ROOTR,ROOTI,IER) 

IF  (IER  .NE.  0)  WRITE  (*,*)  'ERROR  with  POLRT,  IER=',IER,STOP 
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DO  475  1=1, Kd 

write (2, 120)  Roared) 
WRITE(3,120)  Roond) 
S(I)=DCMPLX(RO0TR(I),ROOTI(I)) 
475      CONTINUE 

MAGPOL=0 
DO  485  1=1, Kd 

IF  (CDABS(S(D)    .GE.  l.ODO)  MA(POL=«AGPOL+1 
485      CONTINUE 

WRTTE(*,*)  'I  of  poles  with  magnitude  <=  l',Kd-MAGPOL 
WRITE  (*,*)  'HIT  ANY  KEY  TO  CONTINUE' 
READ  (MOO)  HEADER 

C    Plot  poles 

NOVERLAY=NOVERLAY+l 

CLOSE(2) 

CLOSED) 

CALL  SUBPLT(NOVERLAY) 

J=0 
K=0 

DO  495  1=1, Kd 

IF  (CDABS(S(D)  .LT.  1.0)  THEN 
WRITE  (*,*)  S(I),CDABS(S(I)) 
J=J+1 
K=K+1 
ENDIF 

IF  (J  .EQ.  20)  THEN 
WRITE  (*,*)  'HIT  ANY  KEY  TO  CONnNUE* 
READ  (*,100)  HEADER 
J=0 
ENDIF 
495   CONTINUE 

WRITE  (*,*)  'Poles  with  magnitude  less  than  one  '  ,K 
WRITE  (*,*)  'HIT  ANY  KEY  TO  CONTINUE' 
READ  (MOO)  HEADER 

GOTO  215 

16   STOP 
END 
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APPENDIX    C.    MATRIX    MULTIPLICATION 


SUBROUTINE  MXMUL(A,B,RA,CA,CB,AB) 

INTEGER  RA,CA,CB 

REAL*8  A(70,70),B(70,70),AB(70,70) 

C    Calculates  matrix  multiplication  of  A  x  B=AB,  where 
C    A=RAxCA,B=CAxCB,AB=RAxCB 

DO  30  1=1, RA 

DO  20  J=1,CB 

AB(I,J)=0.0 

DO  10  K=1,CA 

AB(I,J)=AB(I,J)+A(I,K)*B(K,J) 
10   CONTINUE 
20   CONTINUE 
30    CONTINUE 

RETURN 


END 
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APPENDIX    D.    GRAPHICS    ROUTINE 
SUBROUTINE  SUBPLT(NOVERIAY) 


C 

C    MS-FORTRAN  Program  using  "Grafmatic"  Library  Subroutines. 

C    Plots  a  Solid  Line  and  Optional  Overlay  Plot  for  Comparison. 

C    Written  by  M.A.  Morgan  with  Latest  Update  August  1989. 

C 

C    Default  Printer  is  "IBM  Graphics"  (e.g.  Epson,  OkLdata,  IBM) 

C    With  Plot  Rotated  90  degrees  From  the  Vertical.  "GrafPlus.Com" 

C    May  be  Run  to  Rotate  Plot  Upright  on  Paper  and  to  Use  a  Variety 

C    of  Impact  Printers.  "GrafLaser.Com"  May  be  Run  to  Use  a  Laser 

C    Printer.   See  Graf Plus/Laser  Manual  From  Jewell  Technology. 

C 

C 

CHARACTER*1  YN,YN1,DUM,YN2, SYMBOL, BELL, FEED, FFYN 
CHARACTERM  LINE 
CHARACTER*7  SYMB 

CHARACTER*16  LTIT,CTIT,FNAME,TITIfR,TrrLEI 
CHARACTER*64  TITLE, HCOPY 

REAL  Cm  (70)  ,CRTI  (70)  ,NRTR  (70)  ,NRTI  (70) 

INTEGER*2  N,JROW,JCOL,ISYMl,ISYM2,ITYPEl,rrYPE2,NSCRN 
INTEGER*2  CYAN, GREEN, WHITE, YELLOW, RED, BLACK, BLUE, NTWO 
INTEGER*2  JROWl,JROW2,JOXl,JCOL2,CROSS,KdPLT,I 
INTEGER*2  PURPLE, RUST 
EXTERNAL  XFUN,YFUNP,YFUNN 
LINF^' —  ' 
WHTTE=7 
GREEN=10 
CYAN=11 
YELLOW=14 
RED=12 
BLACK=0 
BLUE=1 
NTWO=2 
PURPLE=5 
RUST=6 
BELL=CHAR(7) 
FEED=CHAR(12) 

C  Qear  Screen  and  Put  Up  Introduction  -  on  Blue  Backgound  for  EGA 

C  Only;    Another  Background  Color  is  Possible  by  Changing  "BLUE" 

C  in  the  Calls  to  QPREG  and  QOVSCN. 

CALL  QSMODE(NTWO) 
CALL  OPREG(0,BLUE) 
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CALL  COVSCN(BLUE) 
VRTTE(*,*)  BELL 

NS=1 

NSCRN  =16 
nYPE2=0 

C  Calling  GRAFMATIC  Routines  and  Plotting  Fl  Solid  Line  Graph 

ITYPE1=0 
ISYM1=-1 
NDCTS1=0 
JR0W1=1 
JROV2=350 
JC0L1=  75 
J00L2=  565 
XKIN=-1.2 
XMAX=1.2 
YMD*=-1.2 
YMAX=1.20 
YCVERX=1.115 
XDRG=0.0 
YORG=0.0 
XST=-1.1 
XFIN=1.1 
YST=-1.1 
YFIN=1.1 

25        CALL  QSMODE  (NSCRN) 

CALL  OJ>LOT(JCOIJ,JCOL2,JRCW,JROW2,X^ 
+l,YOVERX,1.5) 
CALL  QSETUP(ND0TS1, CYAN,  ISYM1, RED) 
IF(XFIN-XST  .LE.  9.0)  XMAJOR=0.6 
IF(XFIN-XST  .LE.  6.0)  XMAJOR=0.4 
IF(XFIN-XST  .LE.  3.3)  XMAJOR=0.2 
IF(XFIN-XST  .GE.  9.0)  XMAJOR=(XFIN-XST)/10.0 

MINOR=0 
LABEL=1 
NDEC=2 

CALL  QXAnS(XST,XnN,XMAJOR,MINDRfLABEL,NDEC) 
YMAJOR=XMAJOR 

CALL  QYAnS(YST,YFIN,YMAJCR,MINOR,LABEL,NDEC) 
c         Plot  unit  circle 
A=-1.0 
B=1.0 

CALL  QCURV(XFTJN,YFUNP,A,B) 
CALL  QCURV(XFUN,YFUNN,A,B) 

IF  (NOVERLAY-1  .LT.  1)  THEN 
IF  (NOVERLAY-1  .EQ.  0)  THEN 


142 


WRITE  (*,3)  N3VERLAY-1 

ELSE 

NZERO=0 

WRITE  (*,3)  NZERO 

ENDIF 

ELSEIF  (NDVERLAY-1  .GT.  1)  THEN 

WRITE  (*,3)  NOVERLAY-1 

ELSE 

WRITE  (\4)  NOVERLAY-l 

ENDIF 

3 

FORMAT  (13/  OVERLAYS') 

4 

FORMAT  (13/  OVERLAY  *) 

REWIND (10) 

DO  20  I=1,N0VERLAY 
READ  (10,110)  KdPLT 
READ  (10,100)  TTTLER 
READ  (10,100)  TTTLEI 
OPEN  ( 2 ,  FILE=nTLER) 
OPEN  ( 3 ,  FILE^nTLEI ) 
NKd=KdPLT 
DO  27  J=l, KdPLT 
READ  (2,120)  NRTR(J) 
READ  (3,120)  NRTI(J) 

IF     (DSQRT(NRTR(J)**24NRTI(J)**2)    .GT.  1.1)  THEN 
NKd=NKd-l 
NRTR(J)=0.0 
NRTKJ)=0.0 
ENDIF 
27        CONTINUE 
PURPLE=5 
RUST=6 

WHTTE=7 

GREEN=10 

CYAN=11 

YELLOW=14 

RED=12 

BLUE=1 
IF  (I  .EQ.  1)  THEN 

CALL  OSETUP(NDOTSl, CYAN, ISYM1, RED) 
ELSEIF  (I  .EQ.  2)  THEN 

CALL  QSETUP(ND0TS1, CYAN, ISYM1, GREEN) 
ELSEIF  (I  .EQ.  3)  THEN 

CALL  QSErUP(NDOTSl,CYAN,ISYMl, YELLOW) 
ELSEIF  (I  .EQ.  4)  THEN 

CALL  0SETUP(ND0TS1,CYAN,ISYM1,BLUE) 
ELSEIF  (I  .EQ.  5)  THEN 

CALL  OSETUP(NDOTSl, CYAN,  ISYM1, WHITE) 
ELSEIF  (I  .EQ.  6)  THEN 

CALL  QSETUP(ND0TS1, CYAN, ISYM1, PURPLE) 
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ELSEIF  (I  .EQ.  7)  THEN 

CALL  O^ETUP  (ND0TS1,  CYAN,  ISYtfL,  RUST) 
ELSE 

CALL  QSETUP (NDOTS1, CYAN, ISYM1, RED) 
ENDTJ 

CALL  (7rABL(ITYPEl,KdPLTfNRTR,NRTI) 
20   CONTINUE 

READ  (*,100)  DUM 

GO  TO  40 

HCOPY='HARDOOPY — >  ENTER  P  OR  p' 

CALL  QPTXT(30,HCOPY,RED,25,l) 

CALL  QCM0V(55,1) 

HCOPY=' 

CALL  OJTXT(40,HCOPY,BLACK,25,1) 

IF  (DUM  .NE.    'P*   .AND.  DUM  ,NE.    'p')  GO  TO  40 

CALL  PJ>SCRN 

OPEN  (l,nLE=,PRN') 

WRITE  (1,160)  FEED 
100   FORMAT (A) 
110   FORMAT (12) 
120   FORMAT (E12. 6) 
160   FORMAT ('  \A,\) 
40   CONTINUE 

CALL  QSMODE(NIVO) 

CALL  QPREG(0,BLUE) 

CALL  QOVSCN(BLUE) 

WRTTE(*,*)  NKd, 'points  were  plotted' 

RETURN 

END 

REAL  FUNCTION  XFUN(T) 

XFUN=T 

RETURN 

END 

REAL  FUNCTION  YFUNP(T) 

YFUNP=SCRT(1.0-T*T) 

RETURN 

END 

REAL  FUNCTION  YFUNN(T) 

YFUNN=-SORT(1.0-T*T) 

RETURN 

END 
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